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ABS'ERACT 


In the present work the integral egnation method has 
Been applied to transonic axLsynnnetric flov;. The governing 
differential equation for transonic axisymmetric flow is 
transformed into an integral equation using Green's theorem. 

Using functional relationship for velocity distribution, the 
double integral equation is reduced to an integral equation 
valid on the body surface. The integral equation is then 
converted into a set of algebraic equations. The method of 
parametric differentiation is then applied to these equations 
to get a system of first order nonlinear ordinary differential 
equations. These are solved using Hamming's predictor-corrector 
method upto the bifurcation point . Beyond the bifurcation point 
the steepest descent method is applied to the system of algebraic 
equations. In order to compare our results, the method of local 
linearization has also been studied in the' process. The results 
are computed for a parabolic arc body of different fineness 
ratios at various Mach numbers by the integral equation approach 
and the local linearization method. The results obtained are 
compared with experimaital results and the more accurate 


numerical results 
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CHAPTER 1 


iNTRODTJCTION 


1 o1 GENERAL 

The efforts of aeronautical research are focussed on the 
transonic range hy the practical proLlem of flight in the transonic 
range and of flight vehicles passing through the speed of sound. 

There is a considerable difference in the flow pattern and the forces 
exerted on a moving body at transonic speeds and low subsonic speeds. 

I 

Shock waves are an essential feature of transonic flow. The occurence 
of shock on a moving body leads to a rapid increase in drag coeffici- 
ent with increasing Mach number. 

The term transonic flow designates flows in which the 
velocity in one region is subsonic while in the remaining region it 
is supersonic. The two regions are separated by sonic lines or sonic 
surfaces (see figure 1). In some cases shock wave is formed when the 
velocity exhibits a discontinuity at the boundary between two regions. 
In the lower transonic range when the free stream Mach number is 
slightly less than unity, there are one or more supersonic regions 
embedded in the subsonic flow while in the upper transonic range the 
supersonic flow encloses one or more subsonic flow regions. If the 
free stream Mach number is increased continuously from zero, the 
transonic range begins when the highest local Mach number reaches 
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unity and ends when the lowest local Maeh mmher approaches unity. 

The Prandtl-G-lauert equation for the perturbation velocity 
potential for flow past slender configurations is 

Although equation (l.2,l) gave good results in the study of subsonic 
and supersonic flows it was found to be incapable of treating 
transonic flows. This failure is evidenced by the calculated value 
of growing to such magnitude that it can no longer be regarded 

as a small quantity when compared to Uoa . In order to allow the 
study of transonic flow, terms upto second order were retained in 
the quasilinear partial differential equation for a steady flow of 
a perfect inviscid gas, from which equation (l.2,l) was derived, 
doing so an equation of the form 

(1- ^ <!zz = 1/L( K'+l) bxx (1.2.2) 

is deduced. The above equation in cylindrical cooirdinates can be 
written as 

[(1- + <>rr + ~ '5’r =0 (1.2.5) 

which is the governing equation for axLsymmetric transonic flow. 

The coefficient for in (1.2.3) makes the transonic 
equation both nonlinear and of mixed elliptic -hyperbolic type. If 
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this coefficient is positive, equation ( 1 . 2 . 3 ) is elliptic and 
represents locally subsonic flowi and if it is negative, equation 
(1.2.5) is h-srperbolic and represents locally supersonic flow. The 
mixed character of the flow field may occur with local supersonic 
region embedded in a subsonic flow or with local subsonic regions 
embedded in a supersonic flow. From the basic assumptions of small 
disturbance theoiy it is found that the equation (1.2.5) is valid 
everywhere in the flow field, both fore and aft of shock wave . 

At subsonic and supersonic speeds equation (1.2.3) is of the same 
order of accuracy as equation (l. 2 .l) although more difficult to 
solve because of its nonlinearity. 

Cole and Messiter ( 1957 ) obtained equation (1.2.5) by 
systematic expansion procedure. They presented a system of equations 
for the first, second and higher approximations . By a similar proce- 
dure they obtained the boundary conditions on the axis of the body. 

An expression for the pressure coefficient and drag on the body is 
also derived. 

1 .2 LITERATUBE SUR-yEY 

The term A„ 6 ^ in equation (I.2.3) introduces complica- . 
tions in its solution. Since no general mathematical theory exists for 
such equations* Aerodynamists had to develop their own methods, for 
its solution, Many methods have been used to treat two-dimensional 
transonic flow using equation (l, 2 . 2 ) such as analytical holograph 
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methods, relaxation methods and Oswatitsch's integral equation 
method. For the transonic flow around smooth "bodies of revolution 
the integral equation method would suit very well hut attempts with 
this method are quite disappointing. 

Oswatitsch and Keune (1955) introduced the parabolic 
method. They proposed that the transonic flow equation could be 
approximated for flows with = 1 , past slender bodies of revolu- 
tion, by replacing where K is the acceleration 

parameter, a constant. The resulting equation is linear parabolic 
and has the form of equati-.-n of hoat ccnd.ucti.jji. Accuracy of the 
above method depends on the value of acceleration parameter, K. 

To solve the parabolic equation a characteristic method was developed 
and the value of K which best fits the results of the characteristic 
method is taken. Although the parabolic method is direct and 
simple in application, the results lack accuracy. The piedicted 
flow is in good agreement with the data for fore bodies where the 
assumption is valid. Poor agreement occurs for the aft part where 
the axial velocity gradient is constant and negative. 

Maeder and Thommen (l956) observed that the assumption 
made "by Oswatitsch and Keune (l955) is not close to the actual 
behaviour of the fluid at all speeds than originally suspected. 

They suggested that the value of K should be determined at the point 
of maximm velocity in the incompressible case. The results obtained 
by them are in fair eigireement with the experimental ones. 

,, , ‘ , , , , 
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Heaslet and Spreiter (1956) applied the integral equation 
method to transonic flovir around slender wings and bodies of revolu- 
tion. Special attention was given to conditions resulting faram 
shock waves and to the reduction of relations to special forms nece- 
ssary for discussion of the sonic flow with = 1 . Results obtained 
for cone-cylinders, wings and wing body combinations at a free stream 
Mach number 1 were seen to be in good agreement with experimental 
results. An example has also been given in which the general result 
is applied to calculate pressure and integrated foines on a family of 
slender, elliptic e one -cylinders. 

Speriter and Alksne (l959) introduced the so called local 
linearization method, which grew out of the parabolic method of 
Oswatitsoh and Keune ( '1956) . The idea is to linearize the nonlinear 
transonic flow equation by replacing either <? yy or in the non- 
linear term by a constant X , solving the simplified equation, and 
then introducing different values for X for different points in the 
flow. This procedure can be considered as equivalent in some sense 
to replacing the original nonlinear partial differential equation at 
each point. Results obtained by such a procedure depend on the 
choice of X and must be assembled to determine the final results. 

This step is accomplished by putting the results into such a form that 
a first order nonlinear ordinary differential equation is obtained 
for the streamwise perturbation velocity component ti = , after X 

is replaced by the quantity it originally represented. In many cases, 
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tills equation is of sufficiently simple form that it can he inte- 
grated analytically and the solution expressed in closed form. In 
other oases, the integration must he done numerically, hut the 
equation is of such a form that standard methods can he applied, 
Spreiter and Alksne (l959) applied this method to planar flow past 
thin airfoils and flow past slender bodies of revolution and 
obtained useful results. Of all the approximate methods presently 
available this has been the most versatile and of consistent 
accumcy comparable with theoretical and experimental results. 


Hosokawa (l959» 1962) observed that the parabolic method 
of Oswatitsch and Keune (1956) and later extended by llaeder and 
Thommen (1956) had two drawbacks. Firstly it is seen that the 
calculated flow is not influenced by the mixed flow character and 
the movement of the sonic point with increasing Mach number is 
given incorrectly. He proposed a refinement of the method to 
remove such defects. The basic concept of the method of Hosokawa 
is that the solution for u = is considered to be written in the 
form 

u = ft + ^ ( 1 .2 .4) 

l\j mm 

where is the correction factor and u = Ox satisfies 




( 1 . 2 . 5 ) 


An approximate equation is obtained for ^ with the use of (I. 2 . 4 ) 
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and (1.2,5). Hosokawa suggested that the constant K should he 
determined so that the solution of equation (l .2.5) should he a 
good approximation in the vicinity of the sonic point on the hody 
surface. Thus the sonic point and the value of K is simultaneously 
determined. The occurrence of shock is also predicted and is seen 
that it is closer to the experimental data. Maeder and Thommen 
(1961) pointed out that Hosokawa' s method is not appropriate where 
the flow acceleration at the sonic point is large . A detached 
shook wave appearing upstream of the hody when > 1 , and sonic 
freezing of the flow cannot he derived hy this method, 

Sandeman (1962) applied the local linearization technique 
to solve the choked wind tunnel flow problems. By the uas of local 
linearization method and suitable adjustment 0 f the scale 0 f the 
transonic asymptotic solution used to represent the tunnel wall 
pressure distribution, the correction can he made of the same coder 
of magnitude as that given hy the more reliable results of transonic 
theory over a useful range of p 3 X>file to wind tunnel dimensions. 

The use of the method for axisymmetric flows has added difficulties 
and there appear to he no more reliable theoretical corrections 
available for comparison with the result of these methods. 

Spreiter and Stahara ( 1970 , 1971 ) combined local lineariza- 
tion method with transonic equivalence rule for studying three 
dimensional flows. They developed programs for three dimensional 
axisymmetric transonic flows about bodies of revolution to determine 
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surface and flow field pressure distributions. Programs are also 
developed to calculate the lower critical Mach number and upper 
critical Ifach number. Comparisons of results predicted by the 
above programs show for the axisymmetric bodies studied, the 
predicted surface and flovir field pressure distributions are in 
essential agreement over the fore part of the bodies and diverge 
near the rear part. This discrepancy is attributed to the shock 
wave boundaiy layer interaction and wind tunnel wall interference 
effects. Results also indicate that the order of magnitude of 
disturbances are smaller in the axisymmetric case than the two- 
dimensional case. 

Pink (1971) observed that the local linearization technique 
for transonic flow past slender axisymnEtric bodies gave good results 
at Mach numbers close to one and poor agreement at other transonic 
speeds. He suggested some changes in local linearization technique 
to provide reasonable predictions near the high and low ends of the 
transonic regime while retaining the previous good agreement at near- 
sonic speeds. A revision of Spreiter and Alksne’s (l 959 ) solution 
for accelerating, locally near -sonic flow and an empirical descrip- 
tion of shock wave strength and location are patched to the local 
linearization equations for locally subsonic and supersonic flows. 
Calculated pressure distributions are in good agreement with experi- 
mental data. 

During the past few years some numerical methods have been 
developed to predict inviscid transonic flows, including embedded 
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shocks ahoTit two dimensional airfoils and slender bodies of revolu- 
tion. The most efficient techniques at present use relaxation 
methods to solve finite difference approximations to the governing 
steady state equations. Bailey (1971 ) described a relaxation 
method for the numerical solution of the small disturbance equation 
for flow pa,st a slender body of revolution. His results for sub- 
critical flows agree well in the region near the forward portion 
of the body. The computed values for supercritical flow showed 
slightly more negative Cp as well as a shock location slightly 
upstream of the experimental location. 

Krupp and lihrman ( 1972 ) developed a numerical method to 
solve the transonic flow equation in application to two-dimensional 
and axisymmetric flows. A difference scheme is used to represent the 
derivatives. An analytic expression deri-ved for the far field by 
Guderley (l950) is used as an outer boundary condition. The diffe- 
rence equations are solved by a line relaxation algorithm. Results 
obtained for <1, including cases with shock waves for airfoils 
fl.nd slender bodies of revolution showed good agreement with experi- 
mental data. 

Recently Sedin (l975) described a method for calculating 
axisymmetric sonic flow, which is an extension of the work done by 
Sedin and Bemdt (l970). The idea is to generalise slender body 
behaviour close to idie body to decompose the velocity potential 
into two new functions . A differential system of two coupled 
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equations is then formed. If the equations for the two new functions 

are regarded one at a time, assuming the other function to he known, 

the system can he temporarily interpreted as two paraholio equations. 

The time like variahles for the two equations are then defined in 

two opposite directions. Then the two equations are integrated 

■condition 

alternatively, using the houndary^as initial values for one of the 
functions and the Guderley (1950) far field expansion as initial data 
for the other function. This technique has heen found to he a fast 
and reliable tool for computing axi symmetric sonic flow. The rate 
of convergence is high. The agreement with experimental data and 
other methods of computation is found to he very good. 

1 .5 FRESEMT mm. 

Mich work has heen done in the case of two-dimensional 
flows rather than the axi symmetric case. The integral equation 
method applied to two dimensional flows has proved to give good 
3?e suits. However, the attempts to solve the problem of mixed 
transonic flaw past slender bodies of revolution by the integral 
equation method have heen disappointing. The local linearization 
technique still remains a powerful tool in obtaining approximate 
solution to -the transonic flow equation for the axi symmetric case. 

The local linearization technique as applied today can at best he 
described as an engineering approximation to obtain the solution to 
the transonic flow problems. 
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In oxder to fill the gap an attempt has been made here 
to solve the problem of mixed transonic axisymmetric flow » 
using the integral equation approach. In order to compare our 
results, in the process, the local linearization technique has 
also been studied. In Chapter 2 the governing equations and 
boundary conditions have been deduced with the help of asymptotic 
matching principle. In Chapter 5 the integral equation method has 
been applied to the axisymmetric case and the integral equation 
is converted into a set of algebraic equations. The solution to 
the set of algebraic equations is obtained by the method of para- 
metric differentiation . Chapter 4 deals with the local lineari- 
zation technique for axisymmetric flow problem. Finally in 
Chapter 5 the computational details has been given. The results 
are computed for a parabolic arc body of different fineness ratios 
and at different Mach numbers. The discussion and the conclusions 
drawn are finally presented. 



CHAPTER 2 


PROBLEM FORMULATION 


2 ,1 INTRODUCTION 

The transonic flow over slender "bodies, in general, 
involve curved shocks over tlie body surface. When the upstream 
Mach number is not far from sonic, the vorticity, as well as the 
variation in entropy, produced in the downstream flow due to the 
curved shocks may be ignored for all practical jurposes. Conse- 
quently the entire flow may be considered isoenergetic, homen- 
tropic and irrotational despite the presence of shocks and it is 
then possible to define a velocity potential. 

The small perturbation theory for subsonic and super- 
sonic flow breaks down at transonic speeds. This becomes evident 
from the linearized differential equation for the perturbation 
potential, which in the limit M^ 1 "becomes 


for axisymmetric flow. Thus it vdll not "be possible to satisfy 
the boundary condition of vanishing disturbance at infinity. 
This failure corresponds to a nonuniformity in the expansion 
procedure. It is implicit- in the linearized theory that the 
dimensionless velocity pertur'bations in x and r directions 
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^1x ’'^1r order of mgnitude On the other hand 

an examination of the exact solution of simple supersonic flow 
like expansion around a coirner, [see e.g. Liepmann and Roshko 
( 1956 )] shows that the perturbations in x and r directions (u,v) 
from a uniform sciic stream in the x-direction are related as 

3/2 

V 'V u ' . A similar result holds for v/eak obliq.ue shock waves in 

the transonic range. This sort of order of magnitude relation 

will be required, in the large, for the perturbations of the basic 

transonic approximation. If it is knovTn that a potential exists 

(as is the case) this last requirement has interesting consequence 

that expansion must be carried out in a distorted set of coordinates. 

For example represent the perturbation potential by e(!)-](x,p) where 
1. 

P = e® r, so that the actual velocity perturbations ( scfi ) , 

are of correct order of magnitude. An interpretation of 
this coordinate P is given by saying that in the transonic range, 
perturbations carried by a body in the flow extends much farther 
in the transverse direction (large r) than in the flow direction 
and that a reduced coordinate p must be used to bring these 
disturbances back in via?. All these physical considerations are 
the basis formers general form o f expansion which must be used in 
the transonic case. The expansion procedure for axi symmetric flow 
has been carried out by Cole and Messiter (1957)* 

In section 2 .2 the basic equations and the boundary 
conditions for axisymmetric flow past bodies of revolution 
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have heen discussed. In section 2,3 transonic floxj' equation for 
perturbation potential and consistent houndary conditions have been 
derived by using the asymptotic matching principle [see e.g. Ashley 
and Landahl 1965] • flie expression for the pressure coefficient 
is given and a small disturbance approximation to the shock jump 
conditions have been obtained. In section 2.4 the equations 
obtained in the previous section have been normalized using suitable 
transformations for subsonic Mach numbers. 

2 .2 BA SIC E QU ATIONS MD BOUKDARY CDIIDITIOHS 

We shall consider the flow around a slender nonlifting 
body of revolution defined by 

r = R(x) = T H(x) (2 .1) 


for small values of thickness ratio t , The body axis is aligned 
with the free stream direction. 

For axL symmetric flow, using the conservation equations 
for compressible, invisoid and nonheat conducting potential flow of 
a perfect gas, the differential equation for velocity potential 4 
reads [see e.g. Ashl^ and Landahl ( 1965 )] > 


,2 2 . ,2 2 ,. 

^4xx “ ‘J’r ■'4rr 


+ 


2 

r 4x ” ^^x*^r ^ 


( 2 . 2 ) 




where a is the sonic velocity defined hy 


^^.5) 

Here T^o a^ represent the free stream values of the flow velocity 
an!?. tie velocity of sound. 


Shock waves are a prominent feature of transonic flows; 
Hence additional relations are needed for the transition through 
the shock, d is continuous across the shock wave, while the axial 
and radial component of velocities (u,v) are discontinuous. The 
velocity components across the shock are related by the shook polar 
equation [see e.g. Cole and Messiter (l957)] 


<v2 

V 





Y-1 

- Y+1 


ft?. 


2 

V+1 




--V2 


s.. 


Y+1 1 


(2.4) 


where subscripts 1 and 2 denote the conditions immediately ahead, of 
and behind the shock, a^ is the sonic velocity ahead of the shock 
and is given by 


2 


a 


1 



Y “1 .'^2 
L_ 


'v 2 

+ Yi 



(2.5) 


The boundary condition to be applied at -the body surface 


is one of tangent tlm at the surface. Hence, 
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The other "boundary condition is that all flow disturbances die 
out at upstream infinity. 

The local pressure coefficient can be obtained from 
[see e.g. Liepmann & Eoshko (1965)] 


Y 

Cp = -V ui - - 1 } (2.7) 

^ Yir 2aco 


For small flovir disturbances equation (2.7) can be simplified to 
obtain 


C 


P 







( 2 . 8 ) 


2.5 EERTURMTIOH EQUATIOITS FOR AXISmETRIC FLOY/ 

Y^e seek first an outer expansion of the form 

where e is a small nondimens ional parameter measuring the 
perturbation level and 

p =6r (2.10) 

with 6 being a flxnction of e so that 
6-^0 as e 0 . 

By introducing the expjansion (2.9) into (2.3) and (2.4), we obtain 
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after neglecting higher order terms in e or 6 

o 


I -e(Y+l) ^ ° 


2. 1_ 
P 


Here (l- 0(e) . Thais comparing order of magnitude of each 

term in (2.11), it is noted that non-degenerate equation is obtained 
by setting 


2 • 

5 'b e 


( 2 . 12 ) 


Let us select 6 = e(Y+1 ) • Equation (2.11j then becomes 





(2.13) 

i£- 1 


where K = 2 

e(Y+l) 

(■2.1k) 


K is called the transonic parameter. Equation (2.1^) is basically 
nonlinear and Mach number enters through the transonic parameter K, 

The only boundary condition available for equation (2.13) 
so far is that flow peirturbations must vanish at large distances, i.e, 

dx ’ dp ° “ (2.15) 


The remaining boundary condition belongs to the inner region and is 
to be obtained by matching of the inner and outer flow solutions. 

The inner flow expansion is sought in the form 


+e4!,^(x,r) ] 


( 2 . 16 ) 



Using expansion (2.16) into (2.2), (2.3) and (2.6) we find that 
the inner flow is governed hy 
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i 1 ^ ^ 

~ =0 
^Irr r ^1r 


f ■d\ ^ ® 

V (x,R) = -T" = f— 

1r ^ ’ dx dx 


(2.17) 


with the solution 


X *1 ^ 

e (x,r) = ^ S (x)lnr + g^(x) 


= / t ^ ^ (x)lnr + g^(x)j 


(2.18) 


The inner solution cannot, of course, gi-ge vanishing 
disturbances at infinity since this boundary condition belongs to 
the outer region. V/e therefore need to match the inner and outer 
flow solutions. For this we shall use the asymptotic matching 
principle . 


The inner solution (2 

so that, ^ 


.18) gives j 



2Tr r 


JT 

e 


(^) 


(2.19) 


Since p = 6 r , by application of the limit matching principle to 
the radial velocity component, with the choice, 

2 

e = T 
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we find 

limfp<>-ip] = (^) = ^'(^) (2.20) 

p 0 2ir T 

the boundary condition for the outer flov/. 

The perturbation velocity potential in case of a slender 

2 2/5 

axisymmetriG body is thus of order t , as compared to for 

the two dimensional case [see Ashley and Landahl (l965)]« Since 
the disturbances are of, at least, an order of magnitude smaller 
in the axisymmetric case for a given thickness ratio, one would expect 
the transonic region to be correspondingly smaller than in two- 
dimensional case. 

Introducing the outer expansion (2.9) for in the expression 

3 

for C , Eqn. 2.8 , and neglecting the terms of order £ and higher 

Jr 

we get, 

_ - o 2 o^ 2 2 

Cp--2e(j^^-e " e blx (2.21) 

If, now one defines the perturbation velocity potential 
for outer flow as 

o 

and return to original coordinates (x,r) the folloTang formulation 
of transonic axisymmetric flow problem is obtained. 

[1 ~ - (y+ 1) I (j +-(>_+ = 0 


(2.22) 
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with the boundary conditions 

r ->-o 

- S'(x) 

2 IT 

'J’r' ^ 0 at infinity _ 

The pressure coefficient is given by 

which on the body surface r = R(x) becomes 


(2.25) 


(2.24) 


(2.24a) 

The shook jump conditions can be obtained using the expansions 
assumed for perturbation velocity potential^ [equation (2.9)] in 
equation (2.4). Hence we write the expansions for velocities as 


u . (l) 2 (2) 

~ = 1 + eu^ ^ + e u^ 

00 


V _ 5/2r (1) (2) 

- e [v + « o . ] 


(2.25) 


Substituting (2.25) the shock polar equation (2.4) one obtains 

eV(u^^^- ^ + e(u^^^+ u^”*^) + 0(e^)](2.26) 
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Equating the coefficients of terms of in (2.26), we obtain 

first approximation for the shock polar as 

( 1 )^ 

[U2^''^- [(1- - (T+l)!^^ ^ 

= 0 

Equation (2.27) can -also be written as 
[(1- I^) <u>- — ^ ly? < u^> ]<u> + <V> = 0 

" CO ' ^ d. CO 

where <u >=» U2- 

<V>= T2 -v^ . 

This may be regarded as a snail disturbance approximation to the 
shock o^mp conditions. 

Prom the irrotationality condition, for flow ahead of 
the shock we can write 



which by the use of divergence theorem yields 


(2.27) 


( 2 . 28 ) 


//u^n^dS = // VvjiT.^ dS 

s s 

Similarly for the other side of the shook, we have 


(2.29) 


/Jug n^ dS = /J V2n^ dS 

s s 


(2.50) 
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Substracting (2.30) from (2.29), we obtain 


J/ [< u> 152 -< ■V> n^] dS = 


0 


(2.31) 


Since S is any arbitrary chosen surface, we can write fiom (2.31) 


<u>n. - <Y> n. =0 

li. 1 


i.e, <u>cose 2 '" <Y>cos0^ = 


= 0 


(2.32) 


Eelation (2.32) combined x^rith equation (2.28) gives 


{ ( 1 - 1!^) < u > - < u^ >} cos 6^ +< v> cos 0 2 = 0 (2 .53) 


a second form 0 f the shook jump condition for weak jumps . 


2 .4 EORIIALIZATIOE OF TEMSONIC FL 0 ?f EQUATIONS 

For subsonic free stream Mach numbers, it is convenient 
to introduce the following transformations 


X = X 


? = r(l- Y 

- ^ {y+Yl'^cc 
1 - ]\^ 


(2.54) 


Then the transonic equation (2.22) for axi symmetric flows can be 
expressed as 

. — — -{- ^ -Ip — * ( 2 . 55 ) 

^ XX rr - r ^x xx v 

r 
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The boundary conditions (2.25) then reduce to 


lim (? $-) = 
^■-^0 


(Y+i) s’(x) 


( 1 — 


(2.56) 


^ o_ ^ 0 at infinity 

!X1 3!? 


and the shock condition (2.28) or (2.35) becomes 


- 1-2 - - 2 

[ < u> >] <u> + <V> =0 


(2.37) 


[ < u > - <u >] cose^ + <^> cose, = 0 


(2.38) 


The relation (2.24 3-) for pressure coefficient gives on the body 


surface 


- 2 ( 1 - _ 2 

— i- X " ^ ^ ) 

(y+1) C 


(2.39) 




5.1 HTTRODTJCTION 

In this chapter the goirerning nonlinear partial differential 
equation of transonic flow is transformed into an integral equation 
using Green's theorem for supercritical speeds [M^ 

The integral equation in cylindrical coordinates is simplified using 
a functional relationship for velocity perturbation. This gives an 
integral equation on the body surface. The equation so obtained is 
converted into a set of nonlinear algebraic equations using Gauss 
Qwadrature . 

The solution to these equations can be obtained through 
iterative schemes. Hov^’ever difficulty arises due to the fact that 
it admits multivalued discontinuous solutions. Hence to control 
the development of an iterative solution to the system, the method 
of parametric differentiation is used| the parameter introduced is 
the so called artificial viscosity. The method is similar to the one 
used by Horstrud (1975) "tbe treatment of two-dimensional symmetric 
airfoils . 

5.2 INTEGRAL EQUATIOH lOBMJLATIOF 
APPLICATION OF GREEN'S THEOREM 

We apply Green's theorem to equation (2.55)» which relates 


surface integral to volume integral. The theorem states that the 





following relation, holds "between any two arbitrary functions of 
and ^ having continuous fi rst and second derivatives in the 
volume V bounded by the siirface S 

^//{ L(^) ]dV = //(o||-^||)dS (5.1) 

Sil) 3$ 

where are directional derivatives along the normal n 

drawn inward to the surface S and L is Laplacian operator. 

Let ;}» satisfy 


L( ^Jj) 


ip — - 

XX yy 


4 - ip^^ 

zz 


= 0 


(5.2) 


the well known Laplace equation in three dimensions. The fundamental 
solution of this equation is 


, _ 1 
■ 4 irH 

with E = [ (x - 5 )^ + (y-n)^ + (z- 4 )^] 
In cylindrical coordinates we write, 

- ~ - Sivi © 

y= rcose, z = r 
n = P cos V , ^ = p sin v 


( 5 . 5 ) 

( 3 . 4 ) 


( 3 . 5 ) 


Then ( 3 . 4 ) gives 
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Now let us consider the integral identity (5»1) a^nd apply 
it to the case of axLsymmetric flow. ¥e shall treat the volume 
and surface integrals separately. 
lOJUm IN[PEGRAL 

Prom equation (5*l)j the integral over volume V is 

. ly = /// ['l»L(o) - 0 L(i|-)']dV (5.7) 

V 

Prom opuaticn (3.2) wo have L(<f') = 0 , 

And from aquation (2.35) 

!,( o) = d- i”” s ^ ^) = i ) 

9x 9 X 

where u = 0- . Hence equation (5*7) hocones 

ly = 2 Hh ~ (u ^) dV (5.8) 

V 

We integrate (3 >8) hy parts in’ direction. This 
leads to integral over the hounding surface of the volume V, 
all of which vanish except that over any shock discontinuity. Thus 

Irr ^ U COS0. dS - ~ ///’f’gu^ dV ( 3 . 9 ) 

—2 —2 -*2 

where <u>=U 2 -'a^ is the jump across a shock surface. The 
upstream unit vector normal to the shock surface i s considered as 

n = i cose^ + j cose^ + k cosQ^ 
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and in equation (3-9) 

cose-| dS = Sgn (cos9^)p dv dp = -P dv dp 

Tlras vre obtain for the ■vDiume iiitegral, 

2 ^2 

- ^ ^ ^ -a >p dvl; - y//f o' u dV (3.10) 

V 

THE SdBFACE IHTEGRAL 

The integral over the body saa.rface S includes (l) the 
cylindrical surface at infinity, (2) the surface surrounding 
the singula?: point x,y,z, (5) the surface Sj^ around any shock 
discontinuity, and (4) the body surface. 

(i) Cylindrical surface at infinity 

The integral over the cylindrical surface at infinity can 
be expressed as 



Tliis integral, with f given by equation (5-3) and ’isrith the assuiiq)tion 
that <5 'u for e> 0 at infinity, vanishes, 

(ii) Integral over S surrounding the singular point 

”*** Ir — 

Since ij/ is singular at the field point P(x,y,z) , where 
E. = 0, it must be excluded from the region of integration. If the 
soarface surrounding the field point is taken as the sphere of radius 
R = , then 
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'I' = 


4tt K 


and 


dip 

= 


4^ R 


If w is the solid angle, the elemental area of the surface is 


dS = Rq doj 


Thus integration over the surface Sp then gives 




f 


J ( -2 
o R 


1 )r,2 

R R ^ o 
o o 


= - o(x,r) 


(5.11) 


(iii) Integral over shook discontinuity ^ 

3^ — 

The quantities ^ continuous across any shock 

surface, while the normal velocity is discontinuous there. Thus 
these surfaces are to he excluded from the region of integration 
in order to apply Green' s theorem. The integral over Sp then 
becomes 


CsH -*11)43. -//< II >VdS (3.12) 


where n is the direction of upstream normal to Sp and < — > 

3 ? 34 

denotes the jump in ^ across the shock. Expressing < > 


in the form 
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< _ <0 >00 30-, - «5 > COS0U 

311 5 I p 

= -[ < u> COSG-| +< V > 00302] 

the integral over the shock discontinuity surface , equation (5,12), 
then becomes 


// ( 

% 


an 



) dS 


[< u> cos 6^ +< V > cos^ ] ^ dS (3.13) 


( i'v) Integral over the body surface 

- 2'n’ £ — 


(5.14) 


From Appendix A, we see that in the limit as P'^O, 0 is of the form 


0 = S (x) Inp + g^(x) 


Also from (3 •3) and (3 ‘6) we have 


lim 

P”^ 0 


9p 


r cos( 6-v) 

" P x2 3/ 2 

4ir[(x -5) + r ] 


Hence , 


- 3il; 

lim Op ^ = 0 „ 
p->- 0 


Thus equation (3.I4) reduces to 
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//(- 


9 n 


t a<5 V 
. ^ ) dS = - 




a 5 

lim J j dv d 5 

P -»T 3 o o 


( 5 . 15 ) 


Bow from the houndary condition ( 2 . 36 ) we see that 


lim p 
p 0 


M. 

9P 


(Y+l)liF 

1 - - £ 


s'(c) 

2 ff 


( 5 . 16 ) 


and from ( 3 * 5 ) together with ( 3 * 6 ) we have 


lim xjj = 
p 0 


1 


4ir[(x-5)^ + r^] 


( 5 . 17 ) 


Substituting ( 3 ‘ 1 5 ) s^d ( 3 .I 6 ) in ( 3 .I 4 ) sj^d simplifying we obtain 


^ ( 1 -it) 


1_ f* S'i () 

o/{(x-^f+X^} ^ 

(3«18) 


Finally from eq^uations ( 3 * 11 )> ( 5 - 1 ' 5 ) aiid (5 •'IB} we obtain 


|l - + 


J J[<u> cose. + <v> cos dS 


-(rtl) t 1_/ s'(;) 

( 1 - l£ ) o/{( 5 -E)^+i?} 


( 5 . 19 ) 
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INTEGRAL EQUATION EOS TELOCITY POTENTIAL 


SuLstitating equations (3 *10) and (5*19) ia Green's 
identity (3-l) SS't an integral equation for velocity potential as 

-(y+i) ^ -X, S (5) dC 

5(x,r) = — ./ 


(1-ug) 4^ i/ 5-5)^+^ 


+ // [ (<'a >-<u^>)cos6. +< v> cosft ] dS 

Sp 2 ' 


+ dV 

Y ^2 


( 3 . 20 ) 


Using shock condition (2.58) we see that the integral over the shock 
discontinuity Sp vanishes. Thus equation (5-20) reduces to 

(y+1)1^ . i s'(c) ... _ 2 

$.(x,r) = - ^ j- / ' d? + /// u dV (3.2l) 

1- IC 0 ./ _^2 -2 Y 2 


^ (^-kY 


+r 


which can he written as 


0 (i,^) = L + ///'/'.u dV 

L V ^2 


(3.22) 


(y+1)JI? 


where 0.. = 


I- 1- 

1 I s’(c) 

«L ■ - 3r /■=T=_^ . 


Q ) + r 


(5.25) 


being the linear value of liie velocity potential. 
Jj 
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IWTEGR&L EQUATION K)R PERTURBATION VELOCITY 

The longitudinal perturbation velocity is given by the 
x-T/ise derivative of equation (3.21). After isolating the field 
point by introducing the limits 5 = x + g, we get, 


u(x,r) = 


5 — — 

Oj^(x,r) + lim ™ JJdvpdP [ 
3x e-^ 0 gx 


x-e 


2 35 


dg + / 


2 35 


d5] 


x+e 



+ lim 

£■*• O 


rj/a 


u^(x5P, v) 


e dvP dp 

[ ^ + P^ -2 r 0 o s( 6- v) ] 


"///2 ^ -LOI’) a.V (5.24) 

V 


In the limit as e ^0, the influence function in the 
integrand of the above double integral is effectively a two dimensional 
pulse function at the point P= r , 0 = v and of strength 2 ''*■ [see Heaslet 

(1956)] 


The expression for u then becomes, 

tiCx.r) = u,{x,r) + j(x,r) - /// j dV 

^ .35 

consider the volume integral in equation (3.24) 

•*2 2 / \ 00 00 2.11 *-2 

III } 7,2'*’ ^ ^ dvPdpd5 f ( 5 ,p,v) 


35 ‘ 


-coO 0 

£ 

r 

4: 

35 


(3.25) 


(x-5)^+ "t 4p^-25p cos(6-v) 
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ror axisymmetric flow and small r(r^ o), we obtain 


2 o 2 „ 

J- f \ .00 00 - 

dT.i ;j f(?,p)-i( = 




-)dp d5 


( 3 . 26 ) 


Substituting equation (5.26) in ( 3 . 24 ), we obtain for axisymmetric 
flow 

_2 

u(x,r) = u^(x,r) + ^ ( 2 ;,r) 



CO CO “2 

il 


: \2 2 

5 ) +p 


-) dp ds 


( 3 . 27 ) 


This is then the final form of integral equation for transonic flow 
past slender bodies of revolution. 


5.5 SDIPLIFICATIOI OF IICTEGRAL EQUATION 

In solving the integral equation (5*24) the difficulty 
arises due to the kernel of the double integral. Hence we introduce 
further simplification. We assume a functional relationship 

- - - - - 

u(x,r) = u(x,o) e (3-28) 


between the velocity on the body surface and the velocity directly 
above the surface of the body in a similar mannei as for the two 
dimensional case [see Norstrud (l973)l» This relationship satisfies 
the condition of zero disturbance at infinity. In onier that (5-28) 
also satisfy the boundary condition at the body surface, we determine 
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"b as follows. We differentiate equation (5-28) with respect to r 
and set r = 0, hence 


9 r r=o 


h 


so that h = 


- 'n(x,o) 

9r r=o 


_ •u(x,o) 

(5 — ) 


r=o 


Now (0“-)_ can he determined hy , using the boundary condition 
r=o 

( 2 , 56 ) which gives on the body surface 

(Y+V 


lim 0 — = 

— TT 

r ^o 1 -1 


R (x) 


Thus 


(Y+l)ffi? u(x,o) 

(l-ld) r"(x) 


(3.29) 


Substituting relation ( 5 . 28 ) into ( 5 . 27 ), we get 

- 1 -2 - 

u(x,o) e = u^(x,r) + j u (x,o) e 



- p 

00 .7 

L 


! 

o 


“2P/b 

pe 



. 2 
P 


dP dC 
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•which in the limit as r^-o, gives 

u(x,o) =u^(x,o) +in(J,o)- -i / — E((x-g| ;h) dg (5.30) 

CO 

where E(|x-g| ; h) is the influence coefficient defined as 
(see Appendix B) 

-2p/Td 

e(|:-5| i b) - t i.( f - - ap ) 

95^ ° (x-c)^ + 

= 2tr[{E,'(Z) -IT^’(Z)} +|{H^"(Z) - ll'’(Z)}l 

^ 2!x- ?! 

with Z = — hg — . 

5 4 SYSTM OP ALGEBRAIC EQUATIOITS 

The integral appearing in equation (5*30) is solved hjr 
quadrature, hy apuroximating the range of integration with a finite 
number N discrete intervals and letting the unknown function iT(x,o) 
he represented hy a mean value within each element. For points 
ahead of the leading edge (x <o) or for points aft of the trailing 
edge (x >a) the expression for u(x, o) as described by ( 3 . 50 ) is 
going to be identically zero. Thus it is presumed in effect that the 
influence ^rom points upstream and downstream of the body on the 
value of u(x,o) can be neglected. Thus equation (3.30) can be 
replaced by the following system of nonlinear algebraic equations. 
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1 — 2 X \ . 

'1 0= 1 


-.N 


^i=^. -^2 


-2 
u. . 

0 D 


(5..3i) 


i = 1,2, . , _ H. 


where the inlluence coefficients are integral functions of E 

and are defined as 


ii t. 


/ E d5 

- A5 
^“T 


= ^ E(|x.-x.j;hp 


Rewriting equation (5.3l) obtain 


^i - 


1 - 2 
2 ^i ” 




■il ^ 

3=1 


- 2 

. . u . 
ID D 


= 0 


(3.52) 


1 ■" 1,2,«9oa oR • 


5.5 APRLICATIOE 0 F THE liETHOD OF PARAt>ETRIC DIEEBRERTIATIOR 


The method of parametric differentiation is applied to the 
system of equations (5.3S)» introduce a parameter v , say 
artificial viscosity, such that 


P. 


^(■SiCiil.u) ■ -^L,- M 'ij "j! ' ° 


R 


(3.55) 


X "Iji^****® bU • 



Differentiation of (5 *3 3) is perforned with respect to P , which 
yields a system of first oi^der nonlinear ordinary differential 
equations . 


37 


dp 


- 1 - 2 1 ^ ^ -2 . „ - du . 


(1- pn^) 


[ ^ u. - il ^ + 2 yu. ^ )] (5.54) 


ap 


D=1 

i = 1 ,2, ... .IT. 

Here we have considered u.. independent of P , so that u_. = 

1 

at p = 0 and = ''^j_(p) is unknown. 


i 

1 


The unknown vector u^ can be found out from differential 
equations (3-34) as the solution to the defined Cauchy's problem 
at P = 1 , with the initial vector solution 


u.. = u at p = 0 
1 


(3.35) 


The condition for a unique solution of the system (5.34) 
in the range 0 i P ^ 1 is the non -vanishing of the functional 
determinant or Jacobian of F. [u. (p),p3 Yd.th respect to u. for 
fixed p. The Jacobian matrix of the system is given by 

3(F^ F^) 


J(u^ jP-2> " " ' 'j"^) 


(3.36) 


^(u^ , . . . .^j^) 

The Jacobian matrix (3.36) is positive definite under the condition 
that the elements on the principle diagonal are positi-ve i.e., 


(1- pUj^) > 0 

-1 


u^ < V 




9 
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For a certain value of max [iL ] , however, the Jacobian 

matrix becomes singular and the solution curve of the Cauchy 
problem can split into two or more solutions. Such a point p = p* 
in the interval 0 ^ p ^ 1 for which 

det J(tl^ jU^ , . . o . Ujj) = 0 

can be treated as a bifurcation point and the associated solution 

■Tf 

Ui = u^(y ) will be designated as the bifurcation solution. 

After obtaining numerically a unique and real solution 
u^ = ■u-^(p) to the set of differential equations (3-34)r losing 
Hamming's predictor-corrector method, at p* < 1 , which represents 
within some error bounds, the best possible approximation to the 
condition of mininum- bifurcation, the next step will involve crossing 
of the singular Jacobian matrix. Therefore a continuous dependence 
of the solution on the initial data is no longer possible. 

In order to solve the system (3.55) at p= 1, method of 
steepest descent i s used. The values at the minimum bifurcation 
point is assumed as the initial values for the steepest descent method. 

For subcritical and continuous supercritical flows the 
velocity distribution obtained at the point of minimm bifurcation is 
assumed as the initial velocity distribution. But for discontinuous 
supercritical flow the velocity distribution at the point of minimum 
bifurcation is assumed with a shook discontinuity. This is so because 
we cannot obtain a discontinuous solution by assuming the continuous 
solution as the initial condition for steepest descent method. 



CHAPTER 4 


L OCAL LIMEARIZATIOH METHOD 


4 .1 lETRODUCTIOE 

In this chapter method of local linearization due to 
Spreiter and Alksne (l959) ds.s heen discussed. In section 4*2 
details of the method as applied to axi symmetric flo?/s has heen 
given. Section 4 •3 deals v/ith modification of local linearisation 
technique as suggested hy Pink (1971)- In section 4*4 an empirical 
approximation for shook nave location and its strength is given. 
Section 4*5 gives method of solution of the equations obtained hy 
application of the local linearization technique . 

4 *2 A PPLICATIOl!? OF LOCAL L BISARIZATION TECHhIQLE TO AXISmETRIC 
PLCypS 

The nonlinear differential equation for axisymmetric flow 
is 

'J'x *xx = (4*2.1) 

The local linearization' technique is app3.ied to equation (4.2 .l) for 
locally suhsonic, locally supersonic and near sonic flows. The 
three oases are disoussed separately in the following . 

(i) Locally suhsonic case 

It is convenient to introduce for the coefficient of 
in equation ( 4 . 2 . 1 ), which is then rewritten in the form 
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'E x:c 


+ 


A +6 = 0 

Yjr , Vt-T* 


rr 


(4-2.2) 


vi'here = 1- K' - ku. It is now assaiaed tliat >-, is neithar 

Jii CO jjj 

?,ero nor infinite and that it varies snfficientlj slowly that 
its derivatives can ho d.isregarded so that it can he considered, 
temporarily, as a constant. The prnhlem is now similar to that 
in linearized theory of subsonic flow past slender bodies of revo- 
lution. The solution of equation (4-2.2) subject to appropriate 
B.C. on the body surface and at infinity is thus given by [ see e.g, 
hard, G.N. ( 1955)1 




E 


1_ 

Att 


Z 

/ 


0 


s' (?)dg 

— 2“~ 2 
(x-5) + A^r 


(4.2.3) 


The expression for velocity •Ug follows by differentiation of 
(4-2.5). It can be approximated for points on the surface of 
smooth slender body by 


"E ■ 


lih 

4tt 


In 




1_ 

4'n’ 


/ 




S (g) 


- Cl 




It , 

S (x) 

4TT 


In X- 


E 


+ 


n. 


(4.2.4) 


where subscript i refers to the values for incompressible floT/ 
14, = 0. Differentiation of (4.2.4) yields, 


= llM 

dx 41’’ 


In Xg + 


du^ 


(4.2.5) 
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Uow, 1- - ku is restored in place of so that, in effect, 
the local value for Xg is used at each point and the subscript 
on u is dropped. Then equation (4.2.5) becomes 


du 

dx 


^■21 -ka) 



^(xju) 


(4.2 .6) 


Equation (4.2.6) is a nonlinear ordinary differential equation of 
first order for u on the body surface which is solved using the 
condition 


u = u^ at S (x) = 0 


(4.2.7) 


(ii) Locally supersonic case 

In this case we introduce Xg for the coefficient of 4^ 
so that equation (4.2.1) reduces to 





(4.2 .8) 


where Xg = IT -1 + ku. It is aigain assumed that is neither 
zero nor infinite and that it varies sufficiently slowly that its 
derivatives can be disregarded. The problem is then equivalent to 
that encountered in linearized supersonic slender body theory. The 
solution for d for points on the surface of a smooth slender 
pointed body of revolution is given by [see e .g. WaixL, G'.3S.(1955) ] 


« _ . -l- S (0 cl£ 


(x-r- 1 r 

n 


(4.2.9) 
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The correspondirig expression for Ujj can be approximated for 
points on the surface of a smooth slender body by 


u = gJiEl In 
H 4-t 


^TT ^ 


4 * 


1 . 


Lx / 
0 


X- 5 


dC 


U , 

s (x) 

4 x 


In X- 


■E 


%(x) 


where fg(x) = 


llAinSLA.^ i_ ^ s "(x)- s''(s) 


( 4 . 2 . 10 ) 

(4.2.11) 


fg(x) can be interpreted as the linear value of u for = /2 . 
Differentiation of (4.2.IO) gives 


- ilU). 1 

dx ~ 4 ^ ^ dx 


(4.2.12) 


If, now Mo,- 1+ ku , is restored in placed of Xg so that, in 
effect, the local value for Xg is used at each point and the 
subscript H on u is dipped, equation (4*2.12) becomes 


S ‘ 1“ (>^.-1+ 1“) (4.2.13) 

= F(x,u) 

Equation (4. 2 . 15 ) is a nonlinear first oilier ordinary differential 
equation for u on the body surface which may be solved using 

u = fg at s"(x) = 0 (4.2.14) 


as the starting condition. 
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(iii) Locally nea r sonic case 

The analyses of subsonic and supersonic cases as described 
above have started by introduction of a symbol x foi^ the coeffi- 
cient of rt and the assumption that A is neither zero nor 
infinite and that it varies sufficiently slowly tha,t it can be 
regarded as a constant in the early stages of the analysis. Since 
the results so obtained terminate if x= i-e., if a point is 
approached at \Thich the velocity is sonic, it is evident that 
additional considerations are necessary to study the flows with 
free-stream Mach number near one in which the acceleration from 
subsonic to supersonic velocities is an essential feature. 

In this case is chosen for the coefficient of 
rather than Thus equation (4.2.1) becomes 

jCX. 

dC -1) 0:^. fp (4.2.15) 


where A- = (y +1 


- k 

33; 


It is assum-ed, once again that. 


X- is nonsingular and that it varies sufficiently slowly that it 
can be considered constant in the initial stages of the analyst .. . 


ITow considering Xp positive and constant, v/ith the use of bound.ary 


conditions [equation 


and the Green's 


(2.25) ] 

theorem associated with the L.H.S. of equation (4.2 .15) yields 

“X r^ 

-y. t A-n J- 

4 ir 4, X - 5 -i- o 0 

(4.2.16) 


* - h! ° a? - 1; pv / pie / cp fp 4« 



vrhexe 


e 


(4.2,17) 


p 


4'n'(x -5) 


Now, Up = can be approximated for points on the surface 

of a smooth slender body Ijr 


u_ 


S (x) 1 .. ^p 


S e 


4 TT 


In 


4 tts 


1_ 

4TT 


f S (x)- S (§) 

/ dp 

X - r ^ 

o ^ 


JL J. 

Xp 9x 



(4.2. 18) 


■where C = Euler’ s constant /4( 0.577215. 


If the free stream Mach number is unity fp vanishes 
and Up can be calculated directly. If the free -stream Mach number 
is not unity, the equation (4. 2. 18) is an integral equation. 

However for Mach numbers near unity, it is only necessary to 
approximate locally by substituting X p/k for or au/s^in the 
tripple integral. Poing so and evaluating the integral in (4.2.18), 
we get 




1 - ^ 
(Y+1)1^ 




i^ln-^ 

4ir At! x 


4- 


X I dg (4.2.19) 


If, once again, is restored in place of Xpj the subscript P 

is dropped, and use is made of the following relation between ^ 
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and ^ along the surface of the "body 


d(u). 

dT" 


R 


= . rK.") M 

‘■9X '■ 3r dx 


I n 

r 9 iu , _ s s 

^ 9x ^ 


R 


4‘n’S 


( 4-2 . 20 ) 


a nonlinear ordinary differential equation is olDtained for u on the 
"body surface <, It is 


*1 ” 
u = O + 


(y+i)i 


g (-?-)- In > fs' 

4iT '■ dx 4irS 


] [■ 


(y +1 )l'i!^ooSe^ 

4irx 


1 ) 


, ^ ]. LishAh 

4tt ^ X - C ^ 

o 


Rearranging we obtain 


T n 

du S S 


exp + G- f /2]} 

S 


ix 4.rS 

li -1 

where G(M^,y) = — ^ — 2 
(y +1 )m„ 

Again we see that equation (4 -2 .21) is of the form 


(4.2 .21) 


and fg is as defined in equation (4.2 .11 ) 


du 

dx 


= F(x,u) 


(4.2 .22) 


To obtain a finite welocity gradient at the point where 

tl 

S = 0, it is necessary to assume that 




(4.2.25) 
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As with the cases for subsonic and supersonic flows, the 

value of u at this axial station then serves as the initial 
condition for the solution of equations (4'>2.6) and (4.2.15)" 
At this station, equation ( 4.2.21) becomes 


(”) 

^dx^ '' ^ 
b =0 


S /duN 

S =0 


2^dx ^ - 


S =0 


11 1 


4 ir 


(Y+l)l/^e^ 


ln[- 


(4.2 .24) 


from which the velocity gradient at the starting point is computed. 


The term G arises from a complicated volume integral in 
equation (4. 2. 18). It is evaluated therein by use of approximations 
valid only for Ifech numbers very close to one. Agreement between 
calculated results and expjerimental data is best at Mach number one, 
where G is zero and becomes poor as G increases in magnitude. 

Equations (4.2.6), (4. 2. 15) snd (4.2.21) are first orier 
ordinary differential equations which can be solved using the 
initial ccnditions as specified in the respective sections. The 
coefficient of pressure C on the body surface can be calculated 

.fcr 

by using equation ( 2 .24a) , 

Cp = - 2 u - (g)^ (4.2.25) 

4.5 eim:« s mdbifigation 

Pink suggested that if a revised solution is to be 
meaningful at high transonic speeds G in equation (4.2 .21 ) should 
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i 

approach - » with increasing so that the velocity 

S =o 

given hy equation (4.2.23) would approach tiAot of linearized 

supersonic theory » Similarly, G should approach -(u. ) „ + 

^ S =0 

.4 

^fpr) „ with decreasing transonic Mach numher. The behaviour 
of G with f'o is shown in Pig. 14 . It is seen that very near a 
Mach number of one, G depends on i?® ard and is independent of 
bo6.y shape. Away from one, G depends on body shape, is inversely 
proportional to fineness ratio squared and is independent 
and Y . The simplest function which has this beha.viour is for 

M < 1 , 

00 ' 


f (%“ '• "f 1 - expE -Xo, „ ]} 


S =0 


S =0 

... (4.5.1) 


for M„_ ^ 1 , 


G = - ^(%) „ { 1 - erp[ x=o /( V2 t^) „ ] } ... ( 4 . 5 . 2 ) 


S =0 




S =0 


(l4-l) 

where X„ = = — ^ is the transonic similarly parameter. 

(Y+1)i4t^ 

With G defined by equation (4.5 * 1 ) 02: (4.5.2) the different solutions 
are patched to obtain the final pressure distribution. This has 
been discussed in section 4.5. 


4.4 SHOCK WA-VE APHIOXIMATIOM 

At a Mach number one the shock wave is assumed to be at 

fl 

the after body location ishere S= o . Prom inspection of measured 
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variations of shock wave location with low transonic Ife-ch number 
it was assunBd that shock wave location varied linearly with the 
transonic similarity parameter . It is then interpolated 
linearly from the minimum pressure location at the lower critical 

Tf 

similarily parameter to the point S = o on the after "body at 

The shock wave strength is also empirically determined. 
The shock wave is arbitrarily assumed to produce a downstream 
axial component of Mach number equal to 0.999. Thus the axial 
oomponent of perturbation velocity downstream of the shock wave is 
empirically taken as 


(0.999 ) 



4 .5 METHOD OF SOLUTION 

For purely subsonic flows the differential equation 
(4.2.6) is integrated with the starting condition given by (4«2,7)> 
upto some point near the nose and then returned to the starting 
point Xg and integrated towards the tail. The calculations 
cannot be carried right to the nose or tail as there is a logarith- 
mic singularity at these locations. For purely supersonic flows 

is 

the differential equation (4 .2 .13)/_integrated with the starting 
condition given by (4.2.14). The integration is first done towards 
the nose and then returned to the starting point and continued 


towards the tail. 
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For the cases of purely subsonic and supersonic flows 

the condition at the starting point is sufficient to determine 

a unique solutior. . Bj.t additional considerations are necessary 

for accelerating transonic flows. This is so because the diffe- 

;/ 4 . 2 . 21 ) 

rential equation^ is singular at the point, where S (X) 
vanishes. Consequently, there exist an infinite number of 
integral curves which pass through that point . Of all these 
curves, however, only one is analytic (all derimtives finite) 
and selection of it suffices to determine a unique solution. 

This assures that the solution for vl/T5^ can be expanded in a 

H 

Taylor series in the neighhoarhood of the point where S (x) 
vanishes. The gradient at the neighbouring points of X_ is 
calculated using equation (4.2.21). This is compared with the 
gradient calculated by equation ( 4 . 2 . 6 ). The equation which 
gives smaller velocity gradient is considered and the solution 
is continued to some point near the nose. ITow, the solution 
returns to the point Xg, ?/here the velocity gradient from equation 
(4«2.2 i) is compared with that from equation ( 4.2.1 5)* The 
equation which gi'ves smaller velocity gradient is considered and 
the solution is continued towards the tail. 

Thus the solution starts with transonic accelerating 
flows given by equation ( 4.2.21), which is patched to subsonic 
local linearization solution given by equation ( 4 . 2 . 6 ) towards 
the nose and to supersonic local linearization solution given by 
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equation (4.2.57) towards the tail. The decceleration from 
supersonic to subsonic local linearization may occur either 
through a shock wave as explained in section (4.4) or smooth 
decceleration through sonic velocity whichever is calculated to 
occur fi rst . 

In each of the preceding steps all differential 
equations are integrated using Hamming’s predictor-corrector method, 
with the initial values calculated by Eunge-Kutta method. 



CHAPmi 5 


COMHJTATIO'WAL DETAII.S, RESULTS .AM) PISCUSSIOI? 

5.1 BTTROSUCTI OH 

Tho oompatational details for the method used in solving 
the system of equations (see Cha.pters 5 and 4 ) are gi-ven and 
discussion of results obtained is presented. Results are computed 
for parabolic aro body of revolution of different fineness ratios 
at different free stream Ifach numbers, by integral equation approach 
and the method of local linearization. The results are compared 
with the experimental data given by Taylor (195B) and the more 
accurate numerical results calculated by Bailey (1971) using a 
relaxation procedure. 

5 2 INT EGR AL EQ UA TION APPROACH 
5 “ 2.1 COMHJTATIOHAL DETA ILS 

The flow chart is shown in Eig. 4.Listing of the computer 
program is given in Appendix C, Definition of the variables used 
in the program is given in the listing. The equation for the 
parabolic aro tody of revolution is denoted ly YI^A) . The subroutine 
LINEAR calculates the linear velocity distribution Uffl using equation 
(4.24) with Xj, replaced hjr 1- 1 ^ • applying the Haranlng's 

l.l.T. !'P.NPUR 

nV 


52 


prodic'boi^-corrsc'bor method the initial values are calculated by 
Runge-Kutta method. The initial condition for Runge-Kutta method 
is taken at b == 0 with U = RBR. The subroutine DERIV calculates 
the velocity derivative F. In RERIY subroutine BICOEF is called 
which calculatOvS the influence coefficient Ey given by equation 
(3.27a). In Haaiiiing's predictor-corrector method, first the 
solution vector is predicted which is denoted hy USP. This 
predicted solution is modified and Hamming’s corrector is applied. 

A convergence test is applied at each step so that the values are 
within the prescribed error bound. After calculating the final 
velocity distribution HSN, at different MEW (= v), a check has 

-*'1 

been applied for the Jacobian matrix. If TJSN is less than (AilEW) 
the process of integration is continued, othewirise, the method of 
steepest descent is applied to solve the system. Por continuous 
flows the initial values for the steepest descent method are 
those at which the Jacobian matrix becomes singular. Por dis- 
continuous flows the procedure shown in Fig. 6 is followed. The 
linear velocity distribution and the velocity distribution obtained 
at the bifurcation point is shoifm in the figure. In orier to 
assume a discontinuous solution at the bifurcation point guidance 
from existing results is taken. The assumed velocity distribution, 
which is now discontinuous is also shown in Pig. 6. These values 
are denoted by HHR in the program. The steepest descent method 
is now applied through the subroutine STEEP. In STEEP subroutines 
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FUNVAL and GRAD ane called which, calculate the functional value 
and gradient of functional value of equation (5.35) respectively. 
Prom the bifurcation point the iteration hy steepest descent 
method is carried upto MIMJ = 1 to get the final velocity distri- 
bution. The pressure distribution GPP is calculated at different 
points along the length of the body corresponding to the final 
velocity distribution. 

5 2.2 RESULTS ARP PISOCrSSIOR 

Results are obtained for a parabolic arc body of fineness 
ratio 10 at Ilaoh numbers 0,7s 0*8 which represent subcritical 
oases, and at 0.9 which represents supercritical continuous case. 
Results are also computed for a body of fineness ratio 10 at Mach 
numbers 0.95s 0«975 and 0.990 which represent, supercritical 
discontinuous oases. Results at 1^= 0.8 and 0.95 a body of 
fineness ratio 10, 12, I 4 are also computed. 

The behaviour of influence coefficient E is shown in Pig. 5. 
Pig. 7 shows the comparison of results obtained for a body of 
fineness ratio 10 at Moo= 0.9 with the experimental results of 
Taylor and numerical results of Bailey. It is observed that the 
results obtained by integral equation approach are in good agree- 
ment with the experimental data near the fore and aft of the body. 
But the agreement is not so good in the central region. Pig. 8 
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and Fig, 9 snow inG compaxison of results obtained for a body of 
fineness ratio 10 at Ilach numbers. O.975 and O.99, with the 
expeiimsital and numerical results. It can be seen from figures 
8 and 9 that the results agree well with the experimental 
numerical results for the fore and aft part of the body. The 
location of 1h.e shook is also predicted closely to that given by 
numerical results. The results, however, diverge just ahead of 
the shock . 

This divergence may be attributed to the assume! functional 
relationship [see equation (5-28)] which represents the transverse 
variation of velocity in terms of the velocity on the body surface. 
This approximation although quite adequate for suboritical flows, 
may not be accurate enou^ to give satisfactory results for super- 
critical flovTs when shocks are present [ see ITixon ( 1975)3 • 5?he 
accuracy of the numerical results may be due to the fact that they 
make use of boundary conditions to approximate inviscid flovir in 
open- jet and porous-wall wind-tunnel sections. These take care of 
wall induced interference effects. In the integral equation approach 
we have not considered such effects. 

Figure 10 shows the Mach number effect on a body of fixed 
fineness ratio 10 at Mach numbers 0.7, 0.8 and 0.9, while Fig. 11 
gives the results at Mach numbers 0.95, 0.975 0.99 • It is seen 

that with the increase in llach number the pressure coefficient 
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decreases over the central portion of the "body and the location of 
the shock \Tave moves towards the tail as the Mach numher increases. 

Pig. 12 shows the results obtained at M* = 0.8 for fineness ratios 
10, 12 and 14 showing thickness effect. We observe that as the 
fineness ratio increases, the pressure coefficient decreases over the 
central portion of the body. Prom Pig. 15 it can be seen that as 
fineness ratio decreases the shook location moves towards the tail. 

5.3 LOCiL 1IPEMI2ATI0H METHOD 

5.5.1 . CO imJTATIOML DETAILS 

The flow chart and listing of the program is given in Pig .14a 
and Appendix D respectively. The various variables used are described 
in the listing. To start with we have to solve equation (4,2.24) to 
calculate the velocity gradient at the point X , where S"(x) vanishes. 
This is achieved by method of successive approximations [ see Carnahan 
(1969)] . Wow the velocity at the neighbouring points of X^ is 
predicted using Taylor' s expansion. The solution is continued from these 
neighbouring points to some point near the nose and to some point 
near the tail by integrating differential equations, (4 *2 .6), (4*2.15) 
and (4,2.21) . Por this Hamming's predictor-corrector method is applied 
through the subroutine HMCTHG. In HAlvUlU subroutine SEMI is called, 
which calculates necessary initial values for the predictor-corrector 
method by Runge-Kutta method. In EEMY subroutine MWU is called which 
calculates the derivative P. In HAMUG first the solution vector is 
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predicted v/hioh i'j denoted by yP in the program. YMI represents 
the modifiod value of YP, Hamming's corrector is applied to 
this modified value YIvU . A convergence test is applied at each 
step to keep the errors within a prescribed limit. 

As oxplr-.ined in section 5 of Chapter 4 need to match 
the different solutions for which we compare the gradients. The 
subroutine DHLTi^, calculates the velocity gradients and these 
gradients are compared to decide the solution to be used. The 
variable IICAS indicates the equation being integrated at each step. 
After obtaining the velocity distribution, the pressure distribution 
is calculated. 

5 • 5 *2 BESULTb MTD PISCUSSIOH 

Results are once again obtained for a body of fineness 
ratio 10 at Mach numbers O.7 and 0.8, which are suberitioal cases, 
for 0.9, which is a supercritical continuous case and for 0.95» 
0.975j 0*99 > 1.025 which are supercritical discontinuous 

oases. The results are also computed for purely supersonic Mach 

numbers 1.05, 1*075? 1.10 and 1.20. 

Pigo 14 gives the behaviour of C with Mach number. 

The results of local linearization method for fin^ess ratio 10 
at Maoh mmters 0.7, 0.8, 0.9 ard 0.95, 0.975, 0.990 are compared 

in Fig. 15 ^S. 16 with the experimental data and the results 

obtained by integral equation approach. 
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We see that the results of local linearization is in good agreement 
with the experimental data. Tlie results of integral equation 
approach for suhoritical Mach numhers seem to he in good agreement 
with experimental data, hut the agreement becomes potirr as the 
Mach number increases. The shock location predicted by the integral 
equation approach for M„= O.95 and =.975 is upstream of that 
predicted by the local linearization technique. But forMco = 0,990 
the shook location obtained by the integral equation approach is 
slightly dovmstream of that predicted by the local linearization 
method and the experimental results. The shock location and the 
shock strength in the local linearization method is obtained empe-» 
rically while in the case of integral equation approach the strength 
and the location of the shock is obtained as a solution of the 
problem . 

Pressure distribution for supersonic Mach numbers 1 .01 , 
1.025, 1 .075» 1.10 and 1.20 are also calculated and are presented 
in Pigs. 17 and 18. 

YIe observe that the method of local linearization is more 
versatile in the sense that it can handle the whole range of Mach 
numbers from subsonic to low supersonic and gives reasonably accurate 


results . 
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5 .4 COMCLUEIOHS 

To coiio’'.'udG we see thr.t the integral equation approach 
is 7/011 suited lor solving the prohlem of axisymnietrio transonic flow 
past slender hodies. We have been sucoossful in formulating the 
prohlem in tho same manner as done for two dimensional case. The 
accuracy of the results obtained hy the integral equation approach 
is of the same order as that obtained for the two-dimensional case. 
The location of the shock wave as predicted by the integral equation 
approach is close to the one predicted by the experimental data 
and the more exact numerical results. 


The accuracy of the results obtained by the integral 
equation approach may possibly be improved by improving upon the 
functional relationship betv/een the velocity on the body surface 
and the velocity above the body surface. Again an examination of 
the boundary condition shows that it needs an improvement. The 
exact boundary condition on the body surface can be written as 


lim 
r y 0 




dx 


in the accelerating transonic flow it is seen that 

same ordecr of magnitade as compared to utiity aiid tlie nonlinear 

term i i is retained in the transoric flow equation. Hence 

• . . dH . 

it appears that it is appropriate to retain term like in 

the hoandary condition, and the prohlem should he solTed accordingly. 
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APFEHDIX A 

ASYl-iPTOEIC BEHAVIOUR OF % , AS r->- Q 

The perturhation equation for transonic flow [equation 
(2.22)] can be written as 

(^- 1 ) 

together with 

« |12E1 as r -»'0 (A.2) 

1r r ' 

where F(x) * R(x) r'(x) is the source strength. Prom (a .2), we 
can express the expansion of ^ near r == 0, in the form 

= P(x) In r + g^(x) + ... (-^.5) 

The first two terms of the expansion are the solution to equation 
(A.1) with the R.H.S. equal to zero. P(x) is known but g^(x) is 
an unknoTO function of x. The successive terms in the expansion 
can he obtained hy iteration. Thus using the first term in (A .5) 
i#e* = I*(x) Inr in the R.H.S. of equation (A«l), we obtain 

=*[ ( Y+1) f’(x) In r - Kl/ln r 
1rr r 1r ^ ^ 

= (y+ 1) p'p*'(lnr)^ -KP lux (A 4) 
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This admits solution cf the form 

K) = P(x) In r + g^(x,K) + (r^ log^r) 

=(_[__ p p ) + o(r log r) (a. 5 ) 

Thus 

~ ^ O(r^ln^r) (A, 6) 
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APFEMDE’: B 


SI I'!rLIFICA?IO K OF TIffl ITORAj. TERM IH EQ.TJATIOH r^.27) 


'' -2 


I = /_ E(ix -5h ^)<?C 


(B. 1 ) 


where S(|x - 5 j:,b) 


bS ( / 

BK 0 






(x - 5 )^+ 


ep) 


(B. 2 ) 


Integration in (B. 2 ) can be carried out by using tables of 
integrals [ see. e.g. Gradshteyn and Ryzhik (1965X ]. Thus wo 
obtain 


E(j X - el; b) = b •—[ |(x -5) {h'(Z)- N^(z)} - (x -S) 3 (B. 3 ) 


with Z 


SLszQ. 

b 


Here and N. are the Struve and Neumann functions of first 
1 1 

order respectively. For small Z these functions can be 
represented in series form as [ see e.g. Losch (196O)] 


H,(Z) 


2 

IT 


Z 

2 


24 

2 2 

i .3 .5 


■*’222 
1.5.5- 7 




N (Z) - |Eln(Z/2) + C] ^ _ 

‘ It 1 1 2 » 21 51 


» « 0 



i_ 

2ir 




{1+^ + (l+j + •^)~. 


After carrying out differentiation in (B. 3 ), we obtain 


sCfx-ejj’b) = bE— ii(~ - [ H^(Z)- N^(Z) ] 


+ f(x- 5 )(^f Ch’(z) -n”(z)] 


b/ 1+ i , ^ 4. i 72 O 1 

b{ Y“l 1+1 z ^ + Z - - Z ^ ] 

[ h’(z) - it’(z) ] + 1 z [ 1+ z II + (||)2 ] 


f H/Z) - N^(z)] 


... (B.4) 


Hera, assuming that Z''' 0(l) and “ , < <1, we can neglect 

d5 

the terms containing the dertTratiye of h in eqioation (B. 4 ). Then 
E is obtained as a function of Z = alone. Thus 

E( |J-C1 I b) 5 e(z) - 2irC h|(Z) - h’(Z) J + irZ[ H^(z)- n"(Z) ] 


2tr{[H^’(Z) ^ 1I:|(Z) ] + |C h"(Z) - n’(z). ]> 

... (B.5) 
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- FINAL P^SSSURf: DISTPI3UTION 

• iNFLUSNCii FUNCTinC 
VFLCCITY DFRIVATTV: 

- FINFMSSS RATIU 
STEP Size 

LIMITS OF Ii\:T..GR*TI<';)fl 

“ TOTAL NJMaER OF OIFFuHLNT CASES 
" ySLOCiTY VECTOR 
‘ LINEAR (/'-LOCITY DIS7R,'BUTI0M 

- PReOICTcD SOLUTION 

• FINAL VELrjCITY DISTRieUTION 

VELOCITY DISTRIBUTION ASSUMrC AT BIFURCATION POINT 
' POINTS ALONG THE Ls.N'^G'’’H UF TFE BODY 
MACH NUMBER 


INPUT variables ARL K f YAM , TC , UPR 


25 FORMAT! l-'X,* YAM =*, F6.3 f IL'X » * FINFNFSS RAT 1 0==«=? 13/) 

O;:- FORMAT! 8F?,4) 

33 FORMAT! F5o 3, 12) 

A6 FORMAT! 5F5. 2) 

6‘ FORMAT! IX,I5,7F15.5/ ) 

1-. FORMAT ! 5X f *I 11 )! t =f=X ! I ) * » I IK i =!'Y ! I ) ’i', 1 CXf * YD ! 1 ) * ? ICX, *YOD ! 1) * , 

5RK,=<'U8R!I)’)‘,11X,*R!I )’!‘,11 C »=^CP (I ) ^/) 

131 FORMAT! 2X,16F8. 3) 

FORMAT! IK, !8,3F?.3,4) 

44* FORMAT! 2/., 16F8,3) 

555 F0RMATI2X,13D! 1H<')/) 

556 F0RMAT!2.K,*NY=*,I 3) 

683 FORMAT! SF7o4) 

7 51 FORMAT! 6F7. 3 ) 

2185 F0RMAT!2X,4L =’«',! 3,5 X,=^AM3W =*,F8»4/J 
2? 85 FORMAT! IX, 16F8«4/ ) 

DIMENSION X!25),:PC25),Y!25),YD!25),YC0(25),U8R!25),E!25,25), 


lP{25),U(2i:,2I) fUlJ(25 ),F(25»2lJ jFHHSS )sCPP (2? ) j Z2( 2;:- ) j '^^SL ( 2: j II ) 
:.USP? IS, D ) ,US(2“, 2i) jUSf-jdSjIi )ifS{k5 ,21 ) , UST { 25 , 2 1 ) , FS'H 25, 21 ) , 
UjRR(25) 

C3MMCN/S ‘T77YAM, P A l, TC 

C3M!«!UN/S ! T / UB R , A H i:; W 

C3MM0N/S-T1/:', N,NM 

CaMMCN/Sf-T.-J/NVK 

I^JIIGcR FR 

MKR = Xw 

N=2'. 

N‘4 = N' 1 
C 

C SXPftlSSION FOR PARABOLIC ARC PHuFIL- 

C 

Y3(A} = (TC/2o , 

READ Si:',(F.{! ), 1 = 1, N) 

00 4 II = l,MKR 

lF(MKRoG£o6) R6A3 2:.- , ( .'U I L ) , iL=:. , n! ) 

READ 33, YAM, FR 
TC=l.i /FLOAT {FR} 

PAI=4eD*ATAN do ) 

PRINT 5i"S 
print 25,Y4M,FH 
PRINT SSC 
C 

C CALCULATION OF LI-NLAR VCLOCnY AND PRESSURE DISTRIBUTION 

C 

PRINT 1; 

CALL LINI.':AR( UBR> 

D3 ic- I=1,N 

IF(X{I )orQ«''o) G3 TO i:. 

A=X ( I ) 

15 Y(I ) = (TC/2,)*( !•- {2. *A' l.)=«‘*2) 

YD{n = 2.*TC*(l.-2o»^A} 

YDD( I)=-4®*TC 

16 R( U=ABSCU8R{I >/YDO{ I) ) 

C4={ l.»YAM^YAM)/( 2* 4*YAM=S'YAM ) 

UBR( I) = UBR{n=i‘C4 

CP(I )=-2,’i'UBR( I)-4.*TC#TC=!‘{(1,'-2**X(I ))=«'^2} 

PRINT 6 :,I,X( I ),Y CD , YD{ I ) , YCD (I } ,UBR (I ) ,R( U , CP (I ) 

U3R{ I) = UBR(I )/C4 
i: C3NTINU, 

C 

C INITIAL CONDITIDNS 

C 

69 w H=0.w5 

AHEW=0« C 
f-IVK=i 

00 85 J=2,N 
U( U,1)=UBR{ J } 
use J,l)=U{ J, 1) 



n o o 


Li CONTINUE 
K1 =<% 

K2 = 2;. 

C 

C RJNGi KUTT4 l«!ETH3D 

C 

OQ 110' L = lf4 

C4LL Df=RlV{F,UtA'1£Wf ',1) 

IF(L.fcQ.l) GG TO 91 
DO 92 J=2fN 
92 FS { J,L)=F{ J> L) 

IFCL.EQ.i) GC TO 115 
91 DO 87 J=2sN 
PHI{J)=F( J»L) 

UJ( J)=U{ JfL) 

U{ J»L) = UU{ J) +Te5<' H*F ( J,L ) 

87 CONTINUF 

AMEW=AMcW+H*('..5 

NVK=2 

CALL DERI V{F ,U»A'1F W, S,L) 

DO 9i. J = 2,N 

PHI ( J) = 2.^F{ J, LH-PHI ( J) 

U( Jt L) = UU( J ) +:.o5'i'H=J'F { JjL ) 

9 CONTINUE 

CALL DERIVSF,U,AM;IHi SjL) 

DO 9;T 0 = 2?^ 

PHI<J) = PHl(J) + 2o<‘FCJ»L) 

U( Jf L)=UUt J) +H*F( J,L ) 

CONTIiMUr 

AMFW=AMt;w+H»i‘ ,o5 

CALL DERI V(F,U,AMCW, SfL) 

DO 12U J=2»N 

PHI { J) = (PHI ( J} +F( J,L n/Svi 
U( J* L+1 }=UU( J) +H*PHI ( J) 

US{ JtL+l)=U{ J,L+l } 

USNt J,L+l)=U(J»Ltl) 

12 CONTINUE 

DO 121 J=2tN 
F{ J,L+1) = F( J,L) 

121 CONTINUE 
PRINT 555 
PRINT 21«5iL,AHeH 
PRINT 444j (U (J »L+1 )t J = 2}N) 

11* CONTINUe 

HAMMINGS PRi: Cl COTR'^CGRRECOTR HGTHOD 

015 DO 101 L = KltK2 
AMEW=AMcW+H 
LA=L+l 

IFlL.GE.a) GO TO 122 



OOO l-«iOOO >**000 


OD lv2 '^=2, S'! 

use J»L + i ] = US C J + 2;. *FSt JyL) Fse Jf L' 1 3 +2® *FS ( J, L' 2 ) )/»c 

USP( J,L + 1 )=US( J,L + i) 
lFeL*NE«4) GO TO 1 2 
FS( J,L+1 )=FS (J ,L) 
i 2 C3NT1NU2 
GO TU 

122 DO 1.^2 J = 2fiM 

US ( J f L+ 1 ) =USN ( J? L - 0 ) +4o { 2 . #FS ( J , L ) -F5 { J , L> 1 3 +2 . S U j L- 2 3 3 / 2 . 

i US P(JjL + i3=US{ JjL + i) 

12'^i- NZ=1 
W = 1 

111 IF(L.bQ«^*'3 GC TO U3 

UST ■ MODIFICATION OF PR..OIC'i.-D VALU-'- 
00 104 J=2tN 

124 UST( J,L + 1)=USP{J, L+I 3 + il 2.. /1 21c { US { J ? L 3 • USP { J, L 3 3 
030 IF(NZ«EQ®i3 GO TO ir-C 

LA=L+i 

CALL DER I V ( F S, USf A ME L A 3 
DO lv60 vJ = 2,N 

APPLiCATIUK OF HAMMINGS CQRRVC''a:jR 

use J,L+13 = {9e*US( JjL3“US( J,L- 2 3 + “ IFS { J » L+1 3 +2. *F S( Jj L3“FS{ J» 

5 13 3 3 /8. 

,6 CONTINU 

GO TO l\ 1 
1,.,6 DO 66L J=2,N 
6Ci FSe J,L+13=FS(JtL3 

CALL DERIVeFSjUSTj sAMEWf CfLA 3 
DO 1L9 J=2,N 

use J»L+13=I). 12 5*( 9.*USeJjL3 US(J,L 23 + 2 « F5 1 J» L + 1 ) +2« =^FSe J » L3 
5FSe J,L“13 3 3 

ly9 USL( J»L+13=USe J»L + 1) 

NZ = NZ+1 
GO TO III 
1 S LA=L+1 

CALL OERiVeFSfUS, AMEWtSf LA) 

03 123 J=2,N 

133 use J,L + 13 = r=® 125*i 9o>S'US { J » L 3 'US ( J* L-2 3 +?o FS ( J* L + i )+2o*FS{ JjL) 

5- FSe Jf L' 1) ) ) 
iy7 continue 

IFeNZ«EQ®l) GO TO 114 

APPLICATION OF COUVERGEf^CL TL'ST 

DO UB J=2tN 

138 IF ( ABS e !•*” e USL C vJ» L+1 3 /US ( J I L +1 3 ) 3 « LT« 'yi,f 1 3 GO TO 118 

PRINT 2185,L,AME(e 


o o r5 o o o o o n o c~. a 


: NY = MY+1 

IF(NY.GT,2) CD ID 11 B 
JZ = 1 

CHECK FDR THf JACOBIAN HATRJ. 

DO 13 J=2,N 

':YZ = AMEW=!'USN{J,L) 

iFC^YZ.Gr:. 1 . ) jz=jz+i 

li CONTINUE. 

IF(JZeGi:c2) GO T3 9 
GO TO 114 
9 L1 = L 

AMEW=AMrW-H 
PRINT 2135,L1,AMEW 
369 READ 757, {URR( J)f J=2,N) 

PRINT 555 

PRINT 2785,(URR{J ),J=2,N) 

PRINT 555 
DO 361 J=2,N 
i URR( J) = URR( J )/C4 
DO 775 J=2,N 
USNC J,L1)=URR( j) 

715 CONTINUi: 

/-‘t AMEW=AMt"W+H 

STEfcPEST D^’SCcNT CALCULATION Y.XHNiUUL 

CALL STF;: PCUSNtLl ) 

AMEWI = l./AMt, Vv 

PRINT 2765,(USN{J,Ll+l), J=2sN) 

PRINT 2185, LI, ANEW 
L1=L1+1 

IF(Ll«Lf,X20) GO ro 27 
GO TO 82c 
DO 116 J=2,N 

116 USL{ J,L+i)=US( J,L+1) 

NZ=NZ+1 

GO TO 1G3-: 

CHECK FOR THE JACOBIAN MATRI . 

lit DO 97G J=2,N 

X<=AMEW«USN( J,L) 

IFIXX.Gt.I,) GO TO 9 

CALCULATION OF FINAL VELOCITY DiSTRIBLTIOH USN 

12B USN( J,L + 1}=US{ J,L+l)"-9./121.*{uS( J,L + 1)-USPC J,L+1) 1 
FS{ J,L+1)=FS{J,L> 

971 CONTINUE 
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PRIr;! 

PRINT 21■3;i,L,A^^.r/rf 
P^INT 5S6t<MY 

PRINT , fUSH( Jf L+1 ) ,J=2,n) 

LA=L +1 

CALCULATION OF DERI^ATIV: OF F’ivAL VTLuCnY USN 

CALL DSRIV(FS,US^,AHFV^i2tLA) 

U 3 6 iU J = 2 f N 

ei-‘ FSN{ J»L + 1>=FS{ JtL+l) 

11 CDNTINUt 

CALCULATION OF P^cSSURT DI STP-T EOT lfj:i CPP 

6 2;. D3 3i1) J=2,N 

C%=C 1*"YAM*YAM )/( 2»A#YAM*YAM ) 

USN( Jf21)=USN( J» 21)*C4 

CPP{ J) = - Zo^USNIJ, 21}- la"2<. =>;;( J ) )^=5'2 ) 

PRINT 555 

PRINT 31“ , J,X{ J),U5N{J»2i),CFPU) 

CONTINUf 
A C3NTINU;:: 

STOP 
i MO 


IBFTC SUBl 


THIS SUBROUTINE CALCULATtS THE VSLOCITY DERIVATIVE' WoRolO AME:W 


E -> INFLUi'NCE COEFFICIENT 
F « VELOCITY DgRIVATIVi: 

U - VELOCITY VECTOR 

SUBROUTINE D£RIV( F,Uf AM2 W, E j L ) 

DIMeNSlON U(25,2l),£I25,25)jR(2f- ) -jXI 25)f Fi 25,21 5 

C0MMCN/SET7/YAM,PAI, TC 

COMMON/S^ Tl/X,N,MN 

C3MMON/St:T3/NVK 

C3MW0N/6RP1/R 

CALCULATION OF IMFLUSNCS COEFFICI 'NT 

CALL INCUEFf E,U,L ) 

DO IV J=2,N 
SUM -&80 
DO 20 K = 2,NIM 



noonoon oorioooo#oon 


Kl = K+l 

IF(L®£Q<, JcANCoNVKol-Q.l) OCl T!5 x 
IFCL.NE.l) GO TG 02^? 

IA = 1 
GO lu 
4 I A — L 

SGM=SUM+; .25 =«'!:( J, {UIK, ) + U( K 1 , 1/,. ) )=!«>>t2/4.4, 

5(LI(K,IA )+UlK],IA) )=««( F{K,I A)+F(Kx, I/U) ) 

G3 TO 2- 
'i 1 !A = L 

SJM=SUM+ 6 25*E { Jj K ( OCKj I A) +U{K1 j JA ) )i^*2/4- 
2 CDNTiNUc 

F( J,L} = (l./( l,-A'i:-W^U{ Jf !A) ) )»!«{ ,5*U{ Jf I a )*U { J , . A ) ■ SUM) 
1 CONTINUE 
RETURN 
END 


BFTC SUB2 


HIS SUBROUTINE GALCULATt-S "h:: IMFLU^^Cr, COGFIC! irNT 


» - iNFLUhNCF CO£FFICIi>rr 
SUBROUTirT iMGO?= (I, U,L} 

DI MENSION X( 25 ),E {25y2S) tR(25) jZ2(25) ,UI 25,21) 

C3MM&N/Si:T?/YAM,PAI,TC 

C3MMON/S ■:T1/;<,N,NW 

C3MM0N/GRPi/R 

CA.= { i,“YAr<1*YAM )/{ 2.4-4YAH--5'YAM ) 

YP3=- 4.«TC 
D3 2ia JK=2, N 
DO 210 KK=2,NN 
KI = KK+1 

Z2(KK)=(XCK1 )+X(KK))«Co5 
Z1=22{KK) 

R{KK)=ABS( {U(KK,L )+U{Kl,L) )/(2«4YP3) ) 

R( KK)=R(KK)4C4 
Z=2o*ABSCK( JK)-Z2 (KK) )/R{KK) 

FDHIZ DENOTES FIRST DERIYATIV^: 3F STRIVE FUNCTION*, 
FDYIZ DENOTES FIRST DERIVATIVE JF NFUf'ANN FUNCTION. 
SDYIZ DENOTES S'cDOIJD DERIVATIVE OF N:; I MANN FUNCTION 
SDHIZ DENOTES SECOND DERiVATlVF 3F STRUVE FUNCTIOM* 
C IS EULf-RS CONST Ai^T 

PI»PAI 

DMM4=1o ... 



■-PS = , •i> ' i 

c=^.'. £77215 
1=1 
J=1 
n = i 

DlJHl=Z/:;U 
SUMl=a* . 

w=r?«', 

DJM3 = 1«' /£o : 

F “ ‘w 

'«r* 

,,W -5 ^ 

r I ^ ’’ ® 

S sJ H ^ ^ ‘ 

Y=G.-. 

Q* 0 ® '»# 

P=“{%25*Z 

XK=-v.25*Z 

SLJM4=0»5 

SUM5=£>of? 

DJM&=w.S 
SUM7=t;.n 
DMM3 = 1®’ 

SiJM9=5. 

CONTINU. 

IF(JouQc3. » GO 13 

DMM.:.:- =DMM. -K 1 .5 /FL OAT { J ) + 1. „ / FLOAT ( J - 1 ) ) 

XX=- AX* ( { . 5 ^ 1 )^^ 2 )/ FLUAT { J* U+l ) ) 

p=-'p*{( *z )**?) /ficrn j*{ j+i ) ) 

DUMi=“DU.Ml*( Z**2) /FLUAT{ I*( I +2 ) ) 

DJM'<=”DUM^*{ Z**2} /FLOAT{ I*( 1+2)} 

0UM6=-0UH6*( { 0.5* Z)**2 3/ FLO/ "' ( J*( J 1) ) 

0UM2=DUM1*FL0AT{ J ) 

SUM1=SUMI+0UH2 
FDH1Z=4,L*SUM1/PI 
DUM4=DUH3*FL0AT{ I *a +1) ) 

SU«3=SUH-i+DUH4 
SDH1Z=SUM3*2.0/P1 
DUM7=DUM6*FL0AT(I ) 

DUM5=DUM6 

Sl)M4=SUM^+DUF5 

SUMS=SUH5+DUP7 

SUM6={2«’'' *( ( AL03( .;c5*Z)+C)*SUM£+SUM4 J42« 0/ {Z**2) 3/PI 

SUM7=SUM’H ( DL!«17*3 Mm 3 / PI 

FDY1Z=SUM6-SUM7 

Y=Y-«-XX*, «25*FL0AT{ Cl +l)*{I+2 3 I 

Q=Q+P*FL0AT( J3 

SUM8 = l2.v*C CALOGC t ,a5*2 3+C 3*Y+SUM5/Z+RR3-4«.<,?/( Z**3 3 3 /PI 
JJ = J+l 

DMM4=DMM4 + {1.0/FLOAT(JJ} +1.' /FLOAT(J4-13 3 



W=:W+ e25*,.v«FLOAr { {I +1 )*{ 1+2 ) 

=W/Pi 

SDYlZ=SUH?'^“SUMl,v 

IF ( { FDHl Z • --F ).L=. l'PS« AMD,:, (SDHlZ s Li;<,2 PS® a?>1Do CFDYiZ "'H)eL 
f: AMD* { SDYi z- I J *L= « i;;p s ) go T; s 7 
t;F = FDHiZ 
!.S = SDH1Z 
'■:H=FUY1Z 
tI=SDYlZ 
J=J + I 
1 = 1 + 2 
n=ii+i 
GO TC 3:' 

7 CONTIiMUe 

CC=FDHIZ‘ FOYIZ 
OD=SDHlZ- SDYIZ 
E( JK,KK) = (CC+0.5<'Z’>=DD)*2* 
fe{ JK,KK)=F{ JK,KK) =«' o25/R(KK) 

21 C3NTINUE 
RETURN 
;.MD 
C 

c 

c 

=i='iBFTC SUB2 
C 

c 

C THIS SUBRCUTIMf. APPLIES 7H.. ST' .P-ST CLSCSNT TECHNIQUE 

C 

C 

SUBRCUTINl STE’='P(U,L) 

DIMENSION U{25,21 )sDS{2S>) 

COMMON/ ScT8/UBR,AMuW 
T=C. 1 
N=24 
C 

c calculation of the functional value 

c 

55 CALL FUNVALi S,Uf L } 

SA = S 

c 

C CALCULATION OF 3RAOI?:M OF FUNCTIONAL VALUE 

C 

CALL GRADCUjCSjL) 

JA=;, 

hr, D3 5 J=2iM 

U( J,L) = U{J,L)“T>!'DS(J ) 

5 CONTINUE 

J4 = J^+l 

CALL FUNVALC S,U,L ) 

IFCSA.Lc.S) GO TO 'i3 



= s 

GO ‘i U 1 

f IF( JA.G'"® I ) GO TO t55 
Jk = l: 

03 :!.G J = 2,N 

U{ J,L) = U( J,L ) + T=^‘DS(J) 

■ C3NTINU:' 

H I, 

c accuracy CH£CKIM3 

IF { GO TO 1 

03 37 IA=2fN 

U{ I AiL+l)=U{ IA»L) 

,.7 CONTINUE 
RETURN 
;*ND 


c 

c 

c 

#IBFTC SU84 
C 
c 
c 
c 
c 
c 
c 


this SUbRiJUTINi^ C ALC bLA\TuS "i'H:, FUNCflCNAL VALUu' 


S ~ FUNCTIONAL VALUE 
SUBRUUTirj.. FUNVAL { Sj U»L) 

DI M3NSI UN U ( 25 1 21) f J 3R ( 2 !5 ) t u ( 2G, Z-j ) f R 12^ ) » a ( ) 

COMMON/ SvjTl/>'»N,Niv| 

COMMON/ GRP 1/B 

C3 M MON/ S FT B / UB R » A MOW 

CALCULATION OF I'iFLU frKCU CO'i FF ICJ'.:HT 

CALL INCOE.FC £»U,L) 

S ® 0 » 

D3 lo J=2tN 
Sl=''»i; 

03 5 K=2iNM 

Sl = Sl+E{ JiK) *C U{K»C) +L(KliL) )=!<*2/4. 

S°sI{UU*L)-UBR(J>+AMEW>i‘( ;.oF*(U( J,L)’»*2+G«25#Sl)))=S^*2 

C3MTINUF 

RETURN 

1;NID 


C 

C 

c 

♦IBFIC SUBS 


C 
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'"HIS SUB9UUT INS :ALCULAT,iS 1 H ' GRAD^■,^T OF F-U^JC•■ TONAL V.auS 
S WoRoTc V'-LOCITV' 


DS GRADIENT OF THE FUNCTIU!^/- L VALUf 
SUBRCUTINI GRAD(J,0S,L) 

DiMiNSIOiJ U( 25»2l ) tDS ( 25 ) , YCD{ 2: ) 25s 25 ) ? C 25 ) ,R{25]f UBR( 2 

COMMON/Si-Tl/X, N,NN 
C3MMON/GRP1/R 
C 3 M « 0 N/ S ■ T 3 / UB R , A M- : W 

CALCULATION CF I'^FLUSNCc COEFFICIENT 

CALL INGOLF(^,Uf L ) 

D3 25 K = £,.NN 

K1=K+1 

DS( K)=i:»oO 

D3 U 1=2, N 

Sl=0.fi 

D3 5 J=2,NN 

Jl=J+l 

Sl = Si+E(I, J ,L)+Lt Ji,L) 

- CONTiHU:; 

m= 2<,*(U{ If L)-'JBR{ I )+ AR'iW>!'{-«' ..T*(U{I,L .25^51) ) 

IFd.BQ.K) GO TO 15 

8=0 . 25#AML fc ( I , < ( U { K, L ) +U { K if L ) ) 

G3 TO 2 

15 8=1.+AMEW*{-U{ I,L ) + . 2 K, K U { K ,L )+U{Kl,L) )) 

20 DS{K)=DS{K)+A«=B 
1 C3MTINU.T 
25 CONTIiMU: 

DS{ 24) = DS{23 ) 

RETURf^ 

SHD 


SUB6 


THIS SUBROUTINE CALCULATES TH?: LINEAR VELOCITY DISTRIBUTION 


SUBROUTINE LINEAR tUBRJ 
DIMENSION X( 25 3,JBR{ 25) 
C3MM0N/ SETT/ YAM, PAIfTC 
C3MM0N/SLT1/X,N,NN 
PI=PAI 

C4=(l-“YAM*YAM)/{ 2«4*YAH=«‘YAH ) 



U U o u u 


OS ?. ( A ) = t o »;‘F a I =«'TC’!‘ T C* ( 1 0 */•, +■* .• '1' A ) 

R.t ( A )=2.-*TC=i' {!\--k^A) 

BETA = 1,' YAM*YAM 
DO Iv I=?5'^4 
A= X n ) 

UBR( 1 ) = DS2 ( A )*( A- 0G( BE'^A ) +3» + -.LiiG ? A ) } ) 

UBR( i )=UBR( I )*■ o25/PT 
UBR( I }=UbR( I )/CA 
1 CONTINUE 
RETURN 

’t'iNTHY 


DATA 
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R':AD FRjAYAM 

J'';4P=i 

T = l®t::/FLO«T( FR) 

A=y,.3*{ l.-'SQRT(3. )/3. ) 

/:3=r, 

Z<=X 

6aLC»EXP(., .577215 ) 

CALL DEVIL(G) 

Y( l,l)= 6«^T=«‘T*X*(26--3o>i‘:-,)-G 
YlNI=Y{'i,i) 

H!4=' i, .1 
T3 = T 
YAM=mMN 
hP—i-.' oti'i,' 1 
Z=U. 21132^9 
CC=yo57T215 

AA=- i2.*TC*TC=t' {1. "2. ) 

CD=AL0G{2.4*YAM=^YAM>!'TC>^'Z’H 1. .P(CC)) 

ce=« i 2 «’«'Tc*Tc*a»- So’f'Z) 

BB = A A*A L OG ( T C’i' { 1 . - Z ) ) " C 2 + A A'i' C U 
P=l.t 
1 = 1 

3: Q=AA*ALOG('P)+BB 

IF{ABS{Q< P)-2P ) 1 
2, P=Q 

J = I + 1 
IF (It a:-’ I) 

I t CONTINU.' 

Y{1,1)=:Y{1»1) + HH<'Q 

;' = X+HH 

H=“r::.4Jy5 

U=l. 

P3 T WT 

PRINT 6l, i’FRjAMN 
PRINT 5t/3 

PRINT I, 1,Y(1,1),X 
PRINT 5... 5 
PRINT 50? 

PRINT 111 
PRINT 507 
K= I 

CALL FAULAOtN) 
l(.. ; CBNTINUt 
STOP 
fcND 


♦IBFTC SUBl 

SUBROUTINE FAULAD(N) 
358 F0RMAT{F6.3) 



F-3RMAT(:v',?6(1H#J ) 

F3 R M aI ( «/■ | ,... f I 2 ) 

F0RHAT{F£.4) 

DIMFNSIDN YP(l,22i ,5),F{lj22 ) , (1 ,22<. ) * YMH 1, 2 >■' ,51 t 
:.Y{ l,22r»l,FH{i,223,5) iCP(22M') 

OIMENSICr-i FSUBd, 22's) ,FTRAN{1,22';5,FSLP2R(1,22 ) 

DlMLNSIOi^ YIN( 1,22, ) ,FIN{ 1,22 ) 

DIMF'NSION X1{25),CPL{2S) 

C3MM0N/GRP4/NCAS, NSTAR 
COMMON/ PATFL/Q, AY AM 
C3MMGN/GRP5/KD 
C3MMON/GRP6/YINI 
COMMON/ S£T 1/ T, AMN,U, PI 

C3MMON/SFT2/DS2, 3S3, DS4, OK, APiNj A PM,A P ,APT 

C0MMCN/S£T5/H 

CDMMON/GRPl/K 

COMMON/Si-T6/X,Y 

C3MMON/SFT8/YMI 

COMMON/ SET 9/ I 

C3MMON/FINK/ZX 

HH= i.'ol 

EPS=0.0. 1 

IF(AMN.Gri.'; .95,AND.AHN.L;,;.l. 23) RdD 6,‘':,SHK 

MM = 1 

1 = 1 

CALL REMYCHtN, I ) 

JJ = 1 

L*3 2. V 3 " 5 i*) 

K— 1 

YP ( J , I +„ , K ) = Y{ J , I 3) + ( 4. / .. . ) 2, =i'F t J, 1 ) F{J,I 1)+2.=«'F{J,I-2)S 

IFiloEQoA) 3{J,4)=:'® 

YMI { J,I+2 ,K)=YP{J,I+l,K)+{ll2./i21*)*E(J»I) 

CONTINUE 

J=JJ 

CALL FUMCiFMjJ) 

YMI ( J,I+l,K+l)=( 1./8. )’i‘{9.=5<Y( Jji > 'Y{ J ,I-’*2)+3.=«'H=i'CFN( JidlyK) 
3+2,^F{J,I )' F{ J,I-1)) ) 

IF(ABS( { YMK J»! + l#Kf 1) YMU J , I +1 , K ) ) / YMI {J,I+1,K+1 ) )«L=«EPS) 

563 TO 15' 

IFtK.GE.3) GO TO 3u 
K=K + 1 
GO TO 43 

£( J,I+1) = C 9./121. )=«<{ YMI{ J, I+ljK+l)“YMI(J,I + l,l) ) 

Y( 3,1+1 )=YMI {3,I + 1,K+1)' 3, i + 1) 

33 = 33+1 

IF(3J.LE.N) GO TO- 5& 

GO RO 35 
H=H*G®5 
DO 3 3=1, N 
Y( J,l)=Y{3,l ) 

GO TO 7-;/ 



.i- C3NTiNU::. 

iF{NCAS.Ni:.2 ) GO TU 777 

CALL DELTA{FSUB, = TRAN,FSUF2Rs I ) 

iFCNSTAk.F.Q.2) G3 TO 666 

IF(FSUB(1,I) .Lr.FTRANd, I) ) lvC/:S = l 

NCAS=2 

GO TU 77? 

66c ^F( FSUPSRCl, I) .LE. FTRAhH Ij ! ) ) r;CAS=v 

r:CAS=2 

77? CONTINU' 

ifinhP.gt.i) go rn 372 

iF(AMN.CQ.ANC.X.3T.SHK ) GO "0 . ? ,1 
G3 TO 372 

Y ( 1 , 1 J = ( , . 9 9 9- AM AM N ) / ( 2 * 4=5^ A MN ) 

X=SHK 
H= D • K 0 5 
NAP = 3 
HCAS=1 

GO TO 7: 

672 C3NT1NUF 
X=X+H 

IF{Xo65o: o97132) GO TC 121 
I = I + i 

CALL VE!>U{F»!) 

CP(I)=-2.«Y{ 1, I)/U- 4o*(T=«'*2)*( do 2o^j!)=<'=«'2J 
PRINT 5 ,/.» Ydfi ),CP(I) t'CAS 

IFCXoGE* , o --le.AND.T.LS.l .22) OC' TO 25. 

MM=MH+1 

IF{MMoGTo2) G3 13 12 
H=“’H 

IF(X.GS.‘' .21) GO TQ 12 '■ 

HH= hH 

YC 1, l)=YlMI+HH*a 
...=a.2113249 + HH 
NSTAR=2 

I2'i CONTINUl; 

IFCX.LT.' .97,AN0.MM. cC.2) GO TC 
IFCX.LT.-c.97.AND. MM.GT.2 ) GO TQ 55 
121 CONTINU; 

PRINT 5-/7 
RETURN 
£*4 0 
C 

c 

c 

=^=I8FTC SUB2 

SUBROUTINE REMYCHjN, IJ) 

S?-;.: F0RMAT{5K,lH=«'f 3C5X»F14.?f3X» IH^) »>:Kf 12) 

DIMENSION FSUBClf 22D)fFTRAN{ l»22? )fFStPSRti»22r) 

■ . DIMENSION YA(lt22:.j)f PHIC 1»22C.) »F(l»22i)tYtlt22&)tZCi»22?) 



DlMv !vjSIL:i ^P(22^') 

C3MMCJN/i,.T £>/>;» Y 
C3M^'CN/3..T1/Tt AHM ,U,PI 
C3MM0N/G^P■•'•|./^CAS, NSTAR 
CPd)* 

P^INT 5: ,Y, Y( 1,1 ),CP{1) ,NCAS 

DO 1=1,3 
DO 12b J=1,N 
YACJ,I)=Y{ j,n 
Z( J,1)=Y( J,I ) 

12 C3NTINU; 

CALL VEHU(F,I) 

DO 13d J = 

I3i - PHI ( J,I 5=F{ J,I ) 

DO l^D 'J = 1,N 

143 Y{ J,I) = YA{ J, «5=^F{ja) 

X=X+v,5>»=H 
CALL VENUCF,!) 

D3 150 J=I,N 

I5v PHIC J,I ) = PHI (J,I ) +F{ J, I) ’i'Z. 

DO 160 J = 1,N 

16 Yt J, I )=Yh( J, I) +H>!'F{J, I )«'vo5 
CALL VENU(F,i) 

DO 2i.,L J=1,N 

2 PHI( J,I )=(PHI( J,n+2o*F( J,I) ) 

D3 : ' :L‘ J = i,N 

i.. Y( J,1) = YA{ J, !) +H«'F(J,I) 

.3=X+ oS’i'H 
CALL VHiiU(F,n 
CO it* J = i , i’'j 

h PHI { J,n = {PHl( J,I )+F(J,I))*,166n6 

DO 131. J = i,N 

1C Y(J,!)=YA(J,I)+H*PHI{J,I) 

00 190 J=1,N 
Y( J, I+l ) = Y( J,I ) 

Z( J,i+l)=YCJ,I) 

19. CONTlNUc 

IA=I +1 

CP ( I A ) = - 2o *Y ( 1 , 1 A ) /U'*4c- =S' ! ) =?'{ Ue‘ 2 ) 

P^INT 5;;.D,X, YC 1,1 A5,CP(IA},NCAS 
IF(HCAS.NE.2) GO TO TOO 
CALL DELTA{FSUB,FTRAN,FSUPSR,I ) 
iF(NSTAR,FQ,2) G3 TO 6.'.:. 
IF{FSUB(i,l).LS.FTRAN{l, I)) NCAS=l 
HCAS=2 
GO TO 733 

6(G IF{FSUPfcR{l,I).LE,FTRAN(l,U) DCAS=3 
rCAS=2 

7Ct CONTINUE 
£0 CO NT I NU 3- 

00 222 1=1,4 



o o o 


00 22? J=1#N 

Y{ j,n=i {j,i ) 

22 2 CONTiNU,,:, 

J J = --; 

RETURN 

‘'JD 

C 

C 

c 

#IBFTC SUB^ 

SUBROUTINE V0NU{F,1) 

DIMSivSIUN F( 1, 22D ),y( If 22' ) 

C3MM0N/S::T1/T,AMM,U,PI 

C3 M HON/ S ET 2 / DS 2 , 3 S J» C S ‘.••j DK , A PH j A P M , E P ,^,P T 

C0MMGN/S?T6/)«, Y 

C3 M MON/ GR P4 / NC AS , N ST A R 

C0MMCN/GRP2/0S,DSl,AS2f ASSfRO 

COMMON/ grp:?/ SULC, G 

COMMON/FINK/ZX 

ZK=X 

IF(NCASo£Q,2 ) GO TO 10 
IFCNCAS.FQ.l ) GO TC 2C 
CALL RAJO.VU) 

F< 1, I)=U*{ {US3*APM^0s,2f/Pn + DS:*">.5*AFT/PI--{DS2r' e5) / {PI’S'CI®-:-"' 1 ) 
i2'*.^iT*’*2) + r2, nr*^2)=* j 
( j 3 to 3 

1 CALL KUMAP (I ) 

CALL OFYIL(G) 

F{ if I} = {DS 1 -#AS 3 )/ ( 4 .*PI*US) + ( 2 ./ ( 2 . 4 »i'AMM>f‘AMN*RK^^ 5 ULC ) I’N 
r XP( [4.=J=PI )=5‘{ (Yd ,1) +G vo5=5'{ AS2=«'ALCG(T*{ A ) ) / (2o*f'I ) 
512 o=^=T#T=«'/’= 5 '( 2 o-?,a*;. )) )/ AS 2 ) ) 

GO TO 5, 

2 CALL GAM/;(I) 

F) 1 f D = D« ( 08?.*'" =2 5*A7P/69+'H« Tf; ‘ DS:':»/7 ! +OS:i*So 25*APT /F I 
5+DS2*0,25*(l.»2,*X)/ ■ /J ) ) 

5' CONTINU;:-.. 

RETURN 


*IBFTC SUB4 

SUBROUTINE RAJEEVd) 

DIMENSION Y( 1,223) 

C3MMCN/St.Tl/T, AM?4,U, PI 

COMMON/ StT2/ CS2, 3 S3, DS4 , UK, aPN ,APM,AP ,APT 
Cl3MM0N/SirT6/X,|Y 



#000 #ooo #000 


^PN={ } + UK^Y(l , n 1. ) 

^«P.M=ABS C i>PW) 

^>PM=ALOGl;'.P^j ) 

APT =AL0G (!=«'{ l.-X) ) 

AP=':«0 

RETURN 

S;ND 


iBF'X SUBS 

SUB ROUT IN- KUMARt n 
DIMILNSIOW Y< 1,223) 
C0MM0N/GRP2/CS,0S1,AS2,A5N,R;'; 
CDMMON/GRF")/£ULC, 6 
C3MMQN/Si,Tl/T,AMM,U, PI 
C3MM0N/StT6/X, Y 
BS =4 « # P I *T I * { ( X - X =» *2 } 

OS 1 = 8« *p I «T=5' T^i' { Y. - 2 « * A * •■:+ 2o : =^;‘; =^, ' ) 
AS 2 = 8* * P 6« *X + 6i ) 

AS3=a.#PI#T*T=f'( •6.+12.*>:) 
RX=2o*T^(;v X’S'K ) 

RETURN 

rND 


iBf-TC SUBo 

SUBRUUTl.. - GAMAd) 

DIMENSION Y< 1,223) 

COMMCN/Si.Tl/T, AHM,U, PI 

COMMON/5;:T2/CS2,3S3, DS4, DK, APN,APH,AP ,APT 
COMMON/ Si: T6/X, Y 

DS2=8.*PI=5‘{T=«'>!'2)«' ( l,-6,«::.+6.=4'( ) 

0S3 = a.*PI=i‘(T>»»S'2)«' { 6,+12,’')=X) 

DS4=96o*Pl*( T=S'«2) 

DK = 2*4*iAMN=5‘«2)/J 
APN=(l. (AMN^*2)+-OK=4‘Y( 1, I ) ) 

APN=ABS{APH) 

AP M=ALOGUPN) 

APT-ALOG( {;.=*R=2 ) ) ) 

RETURN 
END 


iBFTC SUB7 

SUBROUTINE FUNC(FN,J) 

DIMENSION YMK 1,22.;, 5) , Y( l,22u), FN( 1, 222,5) 
CDMMON/GRPl/K 



n o n 


C 3 M N I'N/ G R P4 / tc S t M G T ^ R 
C3MMON/GRP2/DSjDSl,AS2f AG'. ,R;; 

C3MMQN/GPP.:; /£:ULC, G 
C3MM0N/S:"-T9/ J 
C3MM0N/S"'Tl/T,AH'i ,U, PI 

CD MMON/ S . -T 1 ' /D S2, D S3 , C S^'s- , DK » i-^K y A PM, A FT 

CDMMOhi/SGTP/H 

CDMMCN/S .T6/Z,y 

COMMON/S: T3/yMI 

CQMMON/.PIfJK/ZX 

X=Z + H 

ZK = X 

IF{NCAS.^-a.2 ) GO TO IS 
IF (NCASoGQol ) GO TO 20 
CALL SI DOU( I ) 

FN( 1,I + 1,K)=U=«'( (DS3*APM* •23/PT)+DS3=«'L.3*APT/P! 

5(DS2*U«i) / {PI*{1. '■>) )-2^j-o*( T^^2)+?2o * )^>. ) 

GO TO 3 

1 CALL FISHU) 

CALL DEVIL (G) 

FMIlfl + ljKJsIOSl* AS3 ) / { Gc *PI *DS ) + ( 2o / ( 2o C H =?= 

sexp ( (4*=i'PI )* {( YMI ( 1, I + l, K) +G- AS2=S'AL0G(T*( 1«-X ) ) / { 2. «Fi ) 

512«4T=!'T=«'X*(2.-3.>i'X)) )/A52) ) 

G3 TO 3 

2 CALL NEWTON! 1) 

FN { 1 , 1 +1 5 K ) = L* ( DS B*D . 2 S’S'APM/ Pi +. .75*0 53/ PI +053*' 25*APT/P I 
b+DS2*L.2:*(l.“ 2.*X)/(FI* :*( lo . ) ) } 

: C3NTINU 

RETURN 
i.ND 


*IBFTC SUBS 

SUBROUTINE SIDOUd) 

DI MENSI ON YM J ( 1, 2 2 , 5 ) , Y ( 1 , 2 2. ) 

COMMON/ s;: tp/f 
C3MMON/GRP1/K 

C3MMON/SfiTli:*/DS2, DS3 »DS4,DK, APN, APM,AFT 

C3MMON/St.Tl/T, AMM,U,PI 

C3MMON/SF:T6/Z,Y 

C3MMCN/SeTa/YMI 

;(=Z+H 

DS2=8.*PI*{T**2) ( l»- 6o*X+6«*( .**21 ) 

DS3 = 8e*Pi*IT**2J* (*“'6®+i26*X) 

0S4=96. *P1*{T**2) 

0K=C AMN**2 )*2.4/U 

#4?^=: ( ( AMN**2 )+DK*YMI {1,1 +1 , X }'.■ ) 

APN=ABS{APN) 

APM==ALOG{APN) 

APT=ALOG{T*{ l.-X) > 



n o ^ o r' 


R”TURt'i 

: VJD 


L 


iBFTC SUB 9 

SUBRCUTaHu NSWTBNId) 

DlMCNSIOf’ YMH 5)=Y(1.2£ J 

C3MMON/GRP1/K 

CD M H ON / S S T i /D S2 1 D£3 i C S4 » K f AP Nj A PM j A FT 

CDMMGN/SCTl/T, AMM,U, PI 

COMMON/Sr:Ti/H 

C3MMON/S rT6/2,Y 

C3MM0N/SF:T8/YMI 

X=Z + H 

DS2=8.*PI>}=«T**2)* ( U- 6o*X+6o*(,:Ti‘*2) ) 

DSa = 8 ®*PI=«'d=>* 2 l*(* 60 + 12 . 

DS4=96. 1^*2) 

0K=2*4=«={AMN'S=’^2 )/U 

APN= ( lo- ( AMN »*2 ) + OK*}* YM I ( I , I + .1 s K ) ) 

APN=ABS {APH) 

APT=ALOG{ {T=i‘*2)*{X (::«=J'2 ) ) ) 

A3M=AL0G( aPN ) 

RETURN 
I NO 


* IBFTC SUBli, 

CU 8 R OUT I N i.: 0 hi T A ( F SU B , FT R A N f F S UP R * i ) 

DI HENSIGI-l FSULMIf 22X ) , FTRAN { 1, 22 , FSCPSR ( 1» 22*' ) , Y { 1 ,22 , ) 
COMMON/ S;:T1/T, AMN,U, PI 
C0MMCN/S:*T6/X,Y 
C3MM0N/GRPA/NCAS, NSTAR 

C0MMcN/S!-T 2/DS2,DS2, DS->, OK, APNm PH, aP , APT 
CD MMON/ GR P 2 / CS , DS 1 , A S 2 , A S 3 , 5d, 

COMMuN/GRP3/£ULC, G 

C3MMON/FINK/ZK 

ZK=K 

IF(NCASoNco2 } GO TO ZC-': 
lF(NCAS.^Q.2.AND. NSTAR.EQ.l) GU TO I'-l 
CALL RAJ?tV(I) 

FSUPER( 1, 1 ) = L*(C3S3=^APM*Oo25/PU +0S3«*So5=<‘APT/PI- { 0S2*Do 5 } /( P I* 
5Cl.-X))-24,*{T=f*>«'2 !=«' #2 

11*='' CONTINUE 

CALL KUMAR (I ) 

CALL DSVILCG) 

FTRAKd ,I ) =IDS1=S'AS3) /{4. ^PI«DS )+ { 2. / ? 2.4*^AHN#AHr-S**RX“^£ULC) 

SfeKPC (4.=»=PI )* {( Y{l,I> +G' AS2>i'ALOG(l*!'{ lo»>; )) / ( 2o > > 

E12o*T*T*X=5'{ 2c-*3«»>!'X n )/AS2)) 
iFINCAS®EQ.2«ANO.N3TAR.e<j.2) GO TC* 2C-L 



tDHTlNU 
CALL GA.4::(I) 

FS U B { 1 , ! ) = i J * ( D S'-j * ' . 2 AP H/ P I + * ? 5 * D S 3 /P i + 0 S 
;'■ +0S2*U. 2'::-’’^ { 1 2» ^‘X )/ { ( 1 , .'..))) 

2' CONTIMU 
»^£TURN 
£MD 
C 

c 

c 


*1 BFTC SUBll 

SUBROUTIIX. FiSHd ) 
dimension YMI{1,22.> S) 

01 MANSION Y{1, 22:i5 

COMMON/ GRP 2/ CS , DS 1 f A S2 f A S ') , R ■; 

COMMON/ GRP3/F.ULCf G 

commqn/sI:;t8/ymi 

COMMON/S.ET1/T, AMM j U» PI 
C3MM0N/S*;T5/F 
COMMON/SeT6/Z* Y 
X=2 + H 

DS=4c*PI*T*T*( ) 

OS 1 =S. * P I T# ( X- ? . * M =*;> 2 . ) 

AS 2= 8. *P I -■i'T^T=»= { 1. - 6. =<‘^+6 . ) 

A S 3 = 8 „ * P 1 ^ T * T { “ S O + 1 2 O =i' J ; ) 

Rn=2®*T*rx'-x*x ) 

FIT URN 


!;M0 


. r;5-4tAPT/PJ 


C 

c 

c 

*IBFTC SUB12 

S'JBRCUTINi- D£VIL(G) 
COMMuN/S'!:, Tl/Tj AM'JjUt PI 

commqn/fink/zx 

x=zx 


2E T A = ( A MN* AMN” 1* ) / ( 2 • 4 'DAMN’S' A MH*T *T ) 
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£ C 2 « " ■^ / ) ) ) ) 
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